arXiv: 1503.0022lv2 [cond-mat.mes-hall] 22 May 2015 


Generalized spectral method for near-field optical microscopy 


B.-Y. Jiang/ L. M. Zhang/ A. H. Castro Neto/’^ D. N. Basov/ and M. M. Fogler^ 

^Department of Physics, University of California San Diego, 9500 Gilman Drive, La Jolla, California 92093, USA 
^Department of Physics, Boston University, 590 Commonwealth Avenue, Boston, Massachusetts 02215 
^Centre for Advanced 2D Materials and Graphene Research Centre, 

National University of Singapore, Singapore 117542, Singapore 
(Dated: May 26, 2015) 

Electromagnetic interaction between a sub-wavelength particle (the ‘probe’) and a material sur¬ 
face (the ‘sample’) is studied theoretically. The interaction is shown to be governed by a series 
of resonances corresponding to surface polariton modes localized near the probe. The resonance 
parameters depend on the dielectric function and geometry of the probe, as well as the surface 
reflectivity of the material. Calculation of such resonances is carried out for several types of ax- 
isymmetric probes: spherical, spheroidal, and pear-shaped. For spheroids an efficient numerical 
method is developed, capable of handling cases of large or strongly momentum-dependent surface 
reflectivity. Application of the method to highly resonant materials such as aluminum oxide (by 
itself or covered with graphene) reveals a rich structure of multi-peak spectra and nonmonotonic 
approach curves, i.e., the probe-sample distance dependence. These features also strongly depend 
on the probe shape and optical constants of the model. For less resonant materials such as silicon 
oxide, the dependence is weak, so that the spheroidal model is reliable. The calculations are done 
within the quasistatic approximation with radiative damping included perturbatively. 

PACS numbers: 68.37.Uv, 71.36.-l-c 


I. INTRODUCTION 

The problem of electromagnetic interaction between 
a material surface and a small external particle is fun¬ 
damental to numerous physical phenomena and spec¬ 
troscopic techniques, including surface-enhanced Raman 
scattering, surface fluorescence, adsorbed molecules spec¬ 
troscopy, and near-held microscopy. From the point of 
view of electromagnetic theory, it is a special kind of 
scattering problem where the scatterer resides in a uni¬ 
form half-space, e.g., vacuum, while the effect of the other 
half-space — the sample — is represented by the sur¬ 
face rehectivity ra{q,u}). The reflectivity may depend 
on the in-plane momentum q, frequency w, and polariza¬ 
tion a = P or S. Far-held optics describes the regime 
q < uj/c. Momenta q ^ w/c, which correspond to in¬ 
plane distances Ap much smaller than the diameter c/w 
of Wheeler’s radian sphere,^ are the domain of near-held 
optics. 

This work is motivated by recent advancements of 
the scattering-type near-held optical microscopy^’^ (s- 
SNOM), which has become one of the leading tools 
for measuring optical response of diverse materials on 
spatial scales as short as 5-20 nm. Thanks to techni¬ 
cal improvements and the development of tunable and 
broad-band infrared sources,the s-SNOM has pro¬ 
vided insights into properties of complex oxides,or¬ 
ganic monolayers,^® graphene, and other two-dimensional 
crystals. 

The schematics of an s-SNOM experiment is shown 
in Fig. 1(a). A sharp elongated probe is brought into 
close proximity of a sample and is illuminated by an ex¬ 
ternal electromagnetic wave with electric held 
Its interaction with the probe creates scattered waves 


^iqp+ik z-iuit^ p _ (^x,y), with arbitrary in-plane mo¬ 
mentum q, including large-g evanescent waves, = 
yjipj/cY — q^ ~ iq. Multiple rehections of these waves 
inside the probe-sample nanogap cause small but impor¬ 
tant changes in the total radiating dipole moment 
of the probe. These changes are detected by measuring 
the far-held scattering signal as a function of the probe 
coordinates. This signal is proportional to the probe po¬ 
larizabilities, 

( 1 ) 

which have the dimension of volume. 

The goal of this paper is to study the properties of 
functions y® and y/ For simplicity, we consider only 
axisymmetric probes. We are especially interested in 
probes of large aspect ratio. In the experiment, strongly 
elongated probes are used because of high longitudinal 
polarizability y®, which promotes an efficient coupling 
between evanescent and far-held radiation modes — the 
“antenna” effect — making the detection of the near-held 
component possible. 

We assume that the length of the probe is much smaller 
than the diameter of the radian sphere c/w, so that the 
scattering problem can be treated within the quasistatic 
approximation. The probe shape we examine the most 
is a prolate spheroid. At hrst glance both of these as¬ 
sumptions are unrealistic because actual probes are not 
spheroidal and their length (typically, tens of /rm) can 
often exceed c/w for w in infrared or optical frequency 
domain. Yet this model was previously found to yield 
quantitative agreement with the s-SNOM experimental 
data for many materials. This apparent agreement can be 
expected in cases where the surface reflectivity ra{q,ui) 
of the sample is not too large, and the aspect ratio of 
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FIG. 1. (Color online) (a) Schematics of an s-SNOM experiment in which a polarizable probe is used to examine a sample 
characterized by the surface reflectivity ra{q). The external electric field Eext incident on the system creates evanescent waves 
inside the probe-sample gap. This modifies the dipole moment p of the probe, which is detectable by its far-held radiation, (b) 
The real-space potential distribution for the hrst four eigenmodes of the probe polarizability computed numerically for a 
spheroidal probe of half-length L — 25a. The axes are the x- and «-coordinates in units of a, the curvature radius of the apex 
of the probe. The probe’s location is represented by the uniformly shaded beige area in the upper left corner of each panel. 


the probe does not vary greatly from one experiment to 
the next. Under such conditions the gross features of the 
s-SNOM scattering amplitude should indeed have only 
a modest dependence on the exact shape of the probe 
and other experimental parameters. However, fine details 
of the scattering amplitude are shape-dependent even in 
this case^® and they may be discerned as the instrumen¬ 
tal resolution improves. Furthermore, for samples with 
high reflectivity, even the gross features become sensitive 
to the shape and size of the probe. To demonstrate these 
trends in this paper we study the longitudinal and the 
transverse polarizabilities in great detail. We will ignore 
the S-polarization reflectivity rs(g, w) because for most 
materials it becomes very small at q ^ Lojc. Hence, 
X'^ are functionals of the remaining reflectivity function 
rp{q,ui) and the probe-sample distance ztip. We show 
that such functionals can be quite complicated, especially 
for strongly momentum-dependent reflectivity typical of 
layered and/or ultrathin materials. Therefore, it is good 
to start with a simpler case of a bulk medium with a 
g-independent reflectivity 

/3{uj) = rp{q,uj), (2) 

so that for a fixed Ztip and w, the probe polarizabilities 
are functions of a single parameter /3. 

It should be clarified that while the absolute reflectiv¬ 
ity may not exceed unity for the radiative modes q < Lojc 
because of energy conservation, for the evanescent ones 
q > w/c it is allowed do so. Large /3’s are indicative of 
weakly damped surface modes in a material, e.g., sur¬ 
face phonons in dielectrics or surface plasmons in metals. 
We use the umbrella term “surface polaritons” for all 
such modes. The energy loss due to evanescent modes 
is governed not by |/3| but by lTo.j3 which must be non¬ 
negative at g > uj/c. (To compute the losses Im/3 needs 
to be integrated over g with a weight that depends on 
the probe-sample distance.In the limit of vanishingly 
small dissipation, Im/3(a;) tends to a i5-function peak 
at the mode frequency. In practice, Im/3 (and conse¬ 


quently |/3|) as high as 10-20 is possible for well-ordered 
crystalline solids, e.g., aluminum oxide AI 2 O 3 possess¬ 
ing sharp phonon modes [Fig. 2(a)]. Therefore, a ro¬ 
bust theoretical formalism must be capable of comput¬ 
ing functions x''{P) the entire upper complex half¬ 
plane. To meet this requirement such a formalism must 
correctly reproduce the analytic properties of functions 
X''{I3)- We adopt a version of the generalized spectral 
method (GSM) in which the total field outside the probe 
and sample is decomposed into eigenfunctions of an aux¬ 
iliary homogeneous problem, and the role of eigenvalues 
are played by the reflectivity /3, the so-called /3-method 
in the terminology of Ref. 21. (Similar formalism is also 
known in the theory of conductivity of heterogeneous me¬ 
dia. ^^’^^) We show that for any probe-sample distance 
Ztip > 0 functions x‘"(/3) are meromorphic. In other 
words, they admit the series representations 

°° ni/ 

:^=Tor||, (3) 

k=o P 

where the sequence of poles /3^ has no accumulation 
points, and so, no upper limit. Additionally, we will show 
that if the probe is made of an ideal conductor and no 
other sources of dissipation are present, then the poles 
/3^ > 1 and the residues R’^ > 0 are real. If the dielectric 
constant ctip of the probe is considered fixed, and f3^ 
depend only on the geometric factors: the probe shape, 
size, and its distance Ztip to the surface. All these re¬ 
sults comply with the general theory of the /3-method 
developed in Ref. 21. 

Both /3^ and grow exponentially with k but their 
ratios remain bounded and satisfy the sum rule 


Here Xo = X‘'(.P = 0) is the polarizability of an isolated 
probe, which does not depend on ztip. These properties 



































3 



a; (cm ^) 


(b) Im P (c) Im P 



FIG. 2. (Color online) (a) Near-field reflectivity Piuj) of bulk AI 2 O 3 discussed in Sec. VI. Whenever the condition Re/3(a;) = 
is met, a local maximum appears in Imx"^. The frequencies of three such resonances are indicated by the dashed lines, (b) 
In the complex plane of p, the poles p'j^ lie on the positive real axis, while real materials trace curves in the upper half plane, 
shown in red. (c) A full electrodynamic treatment predicts that the poles shift into the lower half-plane and an additional 
nonanalyticity in the form of a branch cut [1, 00 ) appears. 


ensure convergence of the series (3) at any P ^ On 
the other hand, if a material-specific P{uj) approaches 
any of /3^, a resonant peak in x‘' and ultimately, in the 
near-held signal, would be observed. 

The divergence of x'' at a given pole implies that a 
nonzero dipole, i.e., free oscillations may exist in the ab¬ 
sence of any external held. Physical intuition about this 
regime is aided by the method of images, according to 
which real charges Qi on the probe interact with their vir¬ 
tual images —PQi inside the sample and for /3 > 1 achieve 
a runaway positive feedback. However, one must keep in 
mind that these eigenmodes arise only in the auxiliary 
problem where the sample is substituted by a hctitious 
material of rehectivity P^. The divergence never actually 
happens in real materials due to their inherent dissipa¬ 
tion, which enters in the form of a positive imaginary 
part in P as shown in Fig. 2(a) and 2(b). The resonances 
are further damped due to shifting of the poles P'^ to 
the lower complex half-plane when radiative corrections 
are considered [Fig. 2(c)], as discussed in more detail in 
Sec. VIII. For a generic probe that ends in a rounded 
tip, the amplitude of the eigenmodes is the greatest near 
the tip, as illustrated in Fig. 1(b) for a spheroidal probe. 
Overall, this physical picture of tip-localized eigenmodes 
is an elegant and economical approach to understanding 
the mechanism of probe-sample coupling. 

The main objective of the present work is to elucidate 
the analytical properties of the coefficients and i?^. 
We focus on the practically interesting case where the 
probe length L is much larger than the curvature radius 
a of the probe tip. We show that for such strongly elon¬ 
gated probes three regimes can be distinguished. The 
first is the short-distance limit ^tip a where the be¬ 
havior of /3^ is universal. We show that it can be derived 
from the known exact solutions for spherical particles 
(Sec. II). The second is the long-distance limit, znp ^ L, 


where the probe acts as a point-dipole and the functional 
form of the resonance parameters is again universal. The 
remaining third regime a < Ztip < L is the most nontriv¬ 
ial one where P^ and depend on the probe shape. 

For all the probe geometries we study the poles P^ 
grows exponentially with k, and so for moderate values 
of P it is permissible to truncate the series in Eq. (3) after 
one or a few leading terms. This truncation is effectively 
done in simplified models®’^^^^^ of the probe-sample cou¬ 
pling, see Sec. IX. However, this simplification may lead 
to qualitatively and quantitatively wrong results at small 
Ztip and/or for large p. The latter characterize highly 
polar materials such as Si 02 ^® (a commonly used sub¬ 
strate) and the already mentioned AI 2 O 3 (an important 
reference material of infrared optics). 

Besides addressing analytical properties of the probe 
polarizabilities, we also discuss methods for their numer¬ 
ical computation. For the simplest case of a momentum- 
independent reflectivity, the calculation can be made vir¬ 
tually instantaneous with the help of Eq. (3) once the first 
few P'p. and are computed and stored. For specihc case 
of a spheroidal probe, this calculation can be further ac¬ 
celerated using the spheroidal harmonics basis instead 
of the standard boundary element method (HEM). Since 
the number of relevant poles and residues is relatively 
small, for further convenience, they can be fitted to an¬ 
alytical forms, see an example for L = 25a spheroidal 
probe in Table I. The speed becomes a crucial consid¬ 
eration if the calculations have to be done repeatedly. 
An important example is extracting optical constants of 
the sample from near-held spectroscopy data by curve- 
htting algorithms.^® One may anticipate to hnd a con¬ 
siderable speed-up if this inverse problem were treated 
using the GSM. The acceleration occurs because the un¬ 
known physical parameter P = P{uj) of the sample and 
the geometric parameters P^ and of the probe stand 
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clearly separated. The GSM also applies for momentum- 
dependent rp{q,uj), e.g., for layered samples; however, 
in the current implementation the speed-up compared to 
the BEM is less significant. 

The remainder of the article is organized as follows. 
In Sec. II we analyze the universal aspects of the short- 
and the long-distance regimes. In Sec. Ill the spheroidal 
probe model is considered. The equations for the poles 
and residues are presented and the results of their numer¬ 
ical solution for the case of a g-independent rp are dis¬ 
cussed. In Sec. IV we explore the effects due to a weakly 
g-dependent surface reflectivity. In Sec. VI we apply our 
numerical method to computing the near-held response 
of bulk AI 2 O 3 , a strongly polar material. In Sec. VII we 
perform the calculation for the same AI 2 O 3 substrate but 
covered with graphene, which is a system with a strongly 
g-dependent rehectivity. In Sec. VIII we discuss the ef¬ 
fects of the probe shape and retardation on these cal¬ 
culations. We also do a similar comparison for Si 02 , a 
less polar material. In Sec. IX we discuss prior theoret¬ 
ical work and close with concluding remarks. Technical 
details of the derivations are summarized in Appendix. 
The source code of our computer program is available as 
the online Supplemental Material for this article. 


II. PROBE-SAMPLE INTERACTION IN 
SHORT- AND LONG-DISTANCE LIMITS 

We start with a qualitative analysis of the short- 
distance regime dehned by the condition ztip 'C a. In 
this limit the structure of the localized polariton modes 
can be understood intuitively by analogy^® to electro¬ 
magnetic modes in an open cavity. The probe-sample 
gap can be approximated by a cavity with height z(p) ~ 
2 tip + (p^/2a) gradually increasing as a function of the ra¬ 
dial position p. For simplicity, let us assume the surface 
reflectivity of the probe is equal to unity, as for an ideal 
conductor. To have free oscillations exist in such a cavity 
the surface reflectivity P of the sample must exceed unity, 
compensating for the exponential decay of the evanescent 
waves. The condition of the self-sustained oscillations is 
/3exp(2zA:^(p)z(r)) = 1. Accordingly, the local radial 
momentum g(p) ~ —ik^[p) = log/3/2z(p). Imposing 
the quasiclassical quantization condition dpq[p) = 
7r[fc -|- 0{1)] for mode number k, we obtain 


where 

a = arccosh -p . ( 8 ) 

It is easy to check that Eqs. (5) and ( 6 ) agree in the 
limit of small a. (Dependence of on a is qualitatively 
similar; however, the residues scale as ^ ka^a^ at 
small a, see Appendix B.) It is reasonable to think that 
the behavior of P^{a) at a ^ 1 should be common for 
any shape ending in a rounded tip. As long as the modes 
are localized at p a, they should be affected weakly by 
the rest of the probe. This hypothesis is supported by 
numerical calculations presented later in this article. 

Consider next the long-distance limit Ztip ^ L. In 
this case the probe-sample interaction can be analyzed 
using the multipole expansion. For the lowest resonance 
fc = 0 it is sufficient to include only the dipole term. The 
dipole moment of the probe is given by p'^ = Xo^tot^ 
where + E^^^ is the total field at the probe 

position and E^^^ is the field induced by the image dipole. 
In particular, EP^ = Pp^Az^^ and = /3pll/84p- 
Solving these equation for p’' and casting the result for 
X'' = P'^/dPext ill fiis form (3), we get 


Po 

- 4z0p/x^ , 

Ro 

- 44p : 

(9a) 


- 8z0p/x{i, 

pll 

- 84p ■ 

(9b) 


For the sphere Xo = so that the last pair of equations 
agrees with the exact result ( 6 ) and (7). The k > 0 res¬ 
onances are dominated by higher-order multipoles. The 
principal dependence of the poles and residues on a is 
expected to be the same as for the sphere, i.e., 

P^, ^ , Rl ~ Ag'^ik + l){k + 2)zl^ if ztip » L , 

( 10 ) 

where = 1/2 and g" = 1. The forms for R'j^ are 
verified numerically in a later section. Equations (9a)- 
(10) imply that in the large Ztip limit the sum rule (4) is 
saturated by the k = 0 mode alone. 

The case of a g-dependent reflectivity can be treated 
similarly. Thus, for fc = 0 one finds^® 


X^iuj.ztip) = 


Xo 


1 - Xo9''{^,Ztip) ’ 

00 

g''{uj,ztip) =c'' /'rp(g,w)e"^‘i"^“pg^fig. 


( 11 ) 

( 12 ) 


log/3fc ~ [fc-I-0(1)] , ztip<.a. (5) 

The mode is localized at distances p < y'ZtipO. The 
validity of this qualitative analysis is supported by the 
exact results for spherical particles. For the = jj part, 
the following compact formulas for the poles and residues 
are available^®^^^: 

/3jl(a) = e(2'=+3)“, (6) 

R^lia) = 4(fc -I- l)(fc -I- 2)a^ sinh^ a , (7) 


0 

where c*- = 1 and = 1/2. Note that the integral 
in Fq. ( 12 ) is dominated by the in-plane momenta g ~ 
l/-^tip, which we assume to be well outside the light cone, 
q ^ uj/c. At Ztip > c/w this condition no longer holds 
and one has to include retardation effects, see Sec. VIII. 

In summary, in this section we presented arguments 
that the limiting case formulas (5) and (9a)-(12) apply 
to perfectly conducting probes of arbitrary shapes. For 
the sphere L = a and for probes of modest aspect ratio 
L > a, these formulas match by the order of magnitude at 
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5.4 5,3. 

TABLE I. Coefficients of the nine-pole rational fits log/3fe = I X] RkoT^Z~^ = X] CiZ' j X] diZ'‘ for L = 25a 

i=0 i=0 2=0 2=0 

and 0.003 < .Z < 10, where Z = Ztipja. The fits for the residues apply only to the first eight poles, fc = 0 through 7. The 
remaining residue Rs, is constrained to obey the sum rule (4). 


k 

as 

a4 

aa 

02 

fli 

ao 

bi 

bs 


bi 

bo 

0 

3 

-36.399 

234.56 

-762.76 

1783.1 

-0.015667 

1 

- 10.345 

83.048 

-417.42 

1522.2 

1 

5 

-25.733 

111.01 

-93.002 

290.11 

0 

1 

-3.4964 

27.841 

11.949 

87.231 

2 

7 

-33.029 

157.55 

-118.29 

1391.3 

0 

1 

-0.9961 

-2.6274 

149.14 

253.32 

3 

9 

-36.251 

173.3 

-40.018 

1879.9 

0 

1 

-0.12625 

-8.2396 

185.03 

246.31 

4 

11 

-47.517 

221.75 

-85.286 

2292.3 

-0.00017314 

1 

-0.3866 

-8.4205 

180.3 

237.17 

5 

13 

-45.678 

233.45 

42.551 

3435.5 

-0.000094847 

1 

1.0728 

-20.329 

253.41 

291.59 

6 

15 

-46.254 

254.23 

223.28 

2547 

0 

1 

-0.75722 

7.941 

175.21 

185.07 

7 

17 

-27.808 

235.93 

770.95 

1961.3 

0 

1 

-1.8957 

34.251 

163.71 

122.51 

8 

19 

-65.583 

251.72 

-308.28 

402.67 

0.000032595 

1 

1.6624 

4.104 

-1.1364 

23.24 


k 

C5 

C4 

C3 

C2 

Cl 

Co 

da 

d2 

di 

do 

0 

3.9999 

303.23 

5141.1 

4811.1 

282.17 

1.4941 

1 

0.77084 

0.023552 

0.000027594 

1 

12.001 

916.75 

17089 

33226 

2371.9 

12.255 

1 

1.4472 

0.11961 

0.0003052 

2 

24.001 

1844.7 

35207 

90005 

36.881 

11.908 

1 

1.9006 

0.012624 

0.00025669 

3 

40 

3166.2 

65393 

288144 

224067 

4417.7 

1 

4.1622 

2.6905 

0.09134 

4 

60.028 

4584.9 

87200 

196974 

5632.8 

28.67 

1 

1.5802 

0.066244 

0.00048749 

5 

84.685 

6304.9 

122561 

226585 

4364.5 

4.237 

1 

1.2714 

0.044389 

0.000077993 

6 

116.26 

8216.9 

166927 

316568 

98948 

843.4 

1 

1.4088 

0.44188 

0.0087351 

7 

146.83 

10354 

214160 

350367 

77249 

334.56 

1 

1.1606 

0.28331 

0.0034808 


Ztip ~ a. However, for strongly elongated probes L 3> a 
an additional intermediate regime a <C Ztip <C L exists 
which requires further study. The simplest example of 
such a shape is the prolate spheroid and we discuss it in 
the next section. 


III. INTERMEDIATE DISTANCES. 

SPHEROIDAL PROBE 

Unlike the problem of a sphere, that of a spheroidal 
probe cannot be solved analytically. However, we can 
take advantage of the separation of variables in prolate 
spheroidal coordinates (Fig. 3), which enables a more 
efficient numerical solution.In this coordinate system 
the spheroid is a surface of constant ^ = L/F = The 
focal length F, the major semi-axis L, the minor semi¬ 
axis W, and the curvature radius a of the spheroid apex 
are related by 

F = \/L2 _ IE2 ^ a = W'^lL. (13) 

This implies = [1 — We assume that the 

major axis of the spheroid is along the z-axis. If the 
distance between the spheroid and the sample is Ztipj the 
sample surface is at z = — L — ztip. 

We consider the quasistatic limit where the scalar po¬ 
tential has the harmonic time dependence oc Its 

spatial part must obey the Laplace equation in the do¬ 
main outside both the spheroid and the sample. There¬ 


fore, it can be expanded into spheroidal harmonics, which 
are products of the generalized Legendre functions of 
the first and second kind PJ^{x) and QY^{x). Here 
m = 0, ±1, ±2 ,... is the z-axis angular momentum and 
I must be greater or equal to |to|. As shown in Ap- 


^ _ U - 



FIG. 3. The prolate spheroidal coordinate system. Contours 
of constant ^ ( 77 ) are confocal spheroids (hyperboloids). The 
unit vector r and polar angle 0 of spherical polar coordinates 
and unit vectors p and z of cylindrical polar coordinates are 
also shown for reference. 
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FIG. 4. (Color online) (a) The first four poles j3'^ of the polarizability for perfectly conducting spheroids. The probe-sample 
distance is parameterized by o = arccosh(l -|- Zti-p/a) [Eq. (8)]. The solid lines are for a strongly elongated spheroid L = 25a, 
the dashed lines are for a nearly spherical one L = 1.01a. The external field is in the 2 -direction, = T. (b) The corresponding 
residues divided by . Poles for different shapes converge at small 2 tip, while residues converge at large 2 tip. (c), (d) 
Similar plots for the external field in the x-y plane, 1 ^ = 11. 


pendix A, the expansion coefficients A^i can be related 
to the charge distribution on the spheroid. For exam¬ 
ple, is proportional to the total oscillating charge 
of the spheroid cx For a passive probe, A'^q = 0. 

The I = 1 terms determine the components of the dipole 
moment induced on the probe: 

Pz = - ipy = ■ ( 14 ) 

For each m the set of coefficients A™; satisfies the infinite- 
order system of linear equations 

OO 

(15) 

i'=i 

where A[p and Hui are defined by Eqs. (21) and (28) be¬ 
low. According to Eq. (14), to find p we need to consider 
only m = 0 and m = 1. The requisite coefficients 5"*; on 
the right-hand side of Eq. (15) are given by 

b\ = -^E^, (16) 

5^1 = 2 (Ej, — lEy) , b ^ 1 = -{Ej; + iEy) . (17) 


If the external field Elext = E^x + Eyij + EzZ is uniform, 
all other b'^i vanish. Once we solve the system (15) for 
m = 0, we can find the transverse polarizability from 


X = 


pZ 


4 , A°i 
9 60i ■ 


(18) 


In turn, the solution for m = 1 would give us A^i and 


II ^ Px - iPy ^ 8 3 ^ 

E^-tEy 9 b\ ' 


(19) 


Equation (15) we wish to solve can be cast in a matrix 
form 


(A"* - H) A™ = b”". (20) 

Matrix A*” is diagonal, A"*;// = AJ^Siy, where 

(-ir 4 
* 2/ -I- 1 etip ~ 1 

and etip is again the dielectric constant of the spheroid. 
If the probe is made of an ideal conductor, etip —t 00 , 


Ctip 


Ql 


‘(Co) 


p, 


I (^ 0 ) A^r(^o) 


( 21 ) 
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then Eq. (21) simplifies to 


The elements of matrix H in Eq. (20) are given by 


4 Qr(go) 

21 + 1 P^i^o) ■ 


( 22 ) 


All these A[" are actually positive numbers because the 
factor (—1)”^ is compensated by the same factor in the 
definition of (5("(Co)- The behavior of A[" at large I is 
approximately exponential, as can be deduced from the 
asymptotic formula 


ao^arccoshCo. (23) 

In Sec. VI we also consider the case where etip is a finite 
positive number. In this case the decay of A[" at large I 
is also exponential but with a different factor in front. 

If the dielectric function of the probe is real and neg¬ 
ative, then A[" can be negative, too. It can also be zero, 
which corresponds to a plasmon (or phonon) polariton 
resonance of the probe. The resonances occur at discrete 
etip values 


_ 

^tip,Z 




1 = 1 , 2 ,..., 


(24) 


see, e.g.. Refs. 33 and 34 (where, in fact, a more gen¬ 
eral case of anisotropic etip is treated). Eor each m, the 
sequence e™p ; asymptotically approaches — 1 as 1 —?► oo. 
The smallest, i.e., the most negative value in each se¬ 
quence is the starting one. It can be alternatively written 
as 


Cp.i = 1 - ^ (25) 

where i/ = T for m = 0, iz = || for m = 1, and L'^ are the 
depolarization factors of the spheroid^® 


All 


(Co" - 1 ) 
1 -L^ 



Co + 1 
Co - 1 



2 


(26) 

(27) 


For prolate spheroids, these depolarization factors obey 
the inequalities 0 < L-*- < < 4, and so e^jp < Ctipp- 

For example, if L = 25a, which we use in our calcu¬ 
lations below, then Co = \/25/24, e°;p = —16.9, and 
eJip I = —1.11. If the probe is made out of platinum or 
iridium, which are common materials for AFM tips, its 
dielectric function can indeed be negative. It is in prin¬ 
ciple possible to achieve plasmonic resonances in such 
probes somewhere in the near-infrared or visible spectral 
range. On the other hand, at mid-infrared frequencies, 
for which we do calculations in this paper, the dielec¬ 
tric functions of such metals are in the range of hundreds 
or thousands. Such probes are very far from any of the 
plasmonic resonances and the approximation of the ideal 
conductor, Eq. (22), can be safely used. We do so in the 
remainder of this Section. 


Hw=2tt J rp{q,uj)Ii+i{qF)Ii,^i{qF)e (28) 
0 

where Iiy{z) are the modified Bessel functions of the first 
kind and 


Zp = Ztip + L . (29) 

As mentioned in Sec. I, the reflectivity rp{q,ui) may have 
strong peaks at the dispersion curves uj{q) of the surface 
polaritons of the sample. In practice, rp(q,uj) is always 
finite, so that the integrand in Eq. (28) is well-behaved 
and exponentially decreasing. A fast method of comput¬ 
ing Flip numerically is explained in Supplemental Mate¬ 
rial. In the remainder of this section we will assume that 
rp{q,u}) is g-independent. We will show that the polariz- 
abilies of the spheroidal probe are meromorphic functions 
as stated in Sec. I. We will also present our analytical and 
numerical results concerning the behavior of their poles 
and residues. 

If rp{q,u}) = /3 = const, then matrix H factorizes H = 
P H and Eq. (20) becomes 

(A™-^H) A"* =b™. (30) 

A particular case of this equation for ztip = 0 was previ¬ 
ously derived in Ref. 32. In general, Eq. (30) implies that 
A™ as a function of (3 has poles that are the solutions 
of the eigenvalue problem 

(A“-4"H)ufc = 0. (31) 

The substitution Uj, = (A™)“^/^Vfe transforms it to 
Vfc=/3^Mvfe, M= (A™)-i/2h(A™)-i/^ (32) 


Since all A(" are assumed to be positive, matrix M is 
real and symmetric, and so its eigenvalues (/S^)”^ are real 
and its eigenvectors can be chosen to be orthonormal. 
Assuming v^, form a complete basis, the solution A"* of 
Eq. (30) can be sought as a linear combination of the cor¬ 
responding Ufc. Taking into account Eqs. (18) and(19), 
we arrive at Eq. (3) with 




= Xo l(vfc)ol 


Xo = 


4 w-H ^3 
9 A™ 


(33) 


where, once again, m = 0 for i/ = T, m = 1 for i/ = 
II, and (vfc)^ is the first component of vector v/.- The 
completeness of the basis entails |(vfc)g|^ = 1, leading 
to the sum rule (4). The explicit formulas for Xo that 
follow from Eqs. (22), (33), (26), and (27) are 


X^ = 


Xo = 


If. 

3^0^ 

2L^ 


^0 


Co + 1 

Co - 1 


1 
Co. 

Co + 4 
Co - 1 


V 


47rL-^ 
n -1 


(34a) 


V 


47rLll 


(34b) 




















where V = (47r/3)L^a is the volume of the spheroid. 
These formulas should be familiar from classical electro¬ 
statics or from the theory of light scattering by small par¬ 
ticles.^^ For strongly elongated spheroid L 3> a, — 1) 
they yield 



£3 

ln(4L/a) ’ 
L^a. 


(35a) 

(35b) 


In Sec. I we stated that the sequence /3^ may not have 
accumulation points. For the present case of a spheroidal 
probe this can be proven directly from the properties of 
matrix M. The first step is to show that the matrix 
elements of H obey the asymptotic bound 


InHw < —{I + 1' + 1) arccosh (36) 


at large I and I'. This can be established using the saddle- 
point integration in Eq. (28). Together with Eqs. (22) 
and (23), this bound ensures that at Ztip > 0 the high- 
order matrix elements of M decay exponentially, 


In Mill < —(^ -|- ^^ -f 1) 


arccosh 


^cosha -I- 


^tip \ 

F ) 



Here a [Eq. (8)] parametrizes the probe-sample distance 
Ztip. Hence, the double series M^, = trM^ is con¬ 
vergent. Considering the identity 


^(/3D-"=trM2<oo (37) 

fc =0 


we see that the accumulation points are ruled out. On the 
contrary, tr diverges at Ztip = 0 and one accumulation 
point does exist: /3 = 1. For the sphere this can be found 
directly from Eq. (7) by setting a = 0. 

In the spherical limit oo an analytical solution 

of our equations exists although it is not obvious. We 
deduced the form of this solution from the method of 
images, see Appendix B. At finite we resorted to solv¬ 
ing the problem numerically. As already mentioned, due 
to an exponential growth of /3(( with fc, only a first few 
of such poles are usually needed for evaluating the po¬ 
larizabilities in question x‘^- To compute such and 
the corresponding we used the following procedure. 
Given L/a and a, we would generate an N x N matrix 
made of the first N rows and columns of the full infinite 
matrix M. We would diagonalize this finite-size matrix 
by standard library routines.^® The obtained eigenvalues 
approximate the first N poles /3((. We would gradually 
increase the matrix size until the poles we are interested 
in would show no variation as a function of N within the 
desired accuracy. The larger L/a and the smaller a, the 
higher N was needed. We found this procedure workable 
as long as N did not exceed about 500. As a rule, the 
higher eigenvalues of larger matrices would either fail to 
reach the accuracy or would show an a-dependence in¬ 
consistent with physical principles. This behavior stems 


most likely from roundoff errors. In principle, one can 
combat them by utilizing higher-precision arithmetic but 
we did not pursue this route. For L = 25a the com¬ 
putation of the first nine poles with at least two-digit 
accuracy was possible for a > 0.08, i.e., Ztip > 0.003a. 
The residues were obtained from the eigenvectors of 
the truncated matrix M using Eqs. (33)-(34b). In the 
interval 0 < a < 0.08 we used the linear interpolation 
between /3 (((q; = 0.08) and /3(((a = 0) = 1. 

The results of these calculations are presented in Fig. 4 
for the first four modes, A: = 0 to 3. The solid lines in 
panels (a) and (c) show /3/- and respectively, as a 
function of a. The corresponding quantities for a sphere 
are shown by the dashed lines. The residues R^/aP are 
plotted in panels (b) and (d). The first nine pole-residue 
pairs of the spheroid for z/ = T have also been fitted with 
an error of 5% or smaller to a combination of elementary 
functions in the range 0.003a < Ztip < lOa. The fitting 
formulas and their coefficients are cataloged in Table I. 
The residue behaves differently from the others be¬ 
cause it was constrained to satisfy the sum rule (4). Using 
these formulas one can find the response with negli¬ 
gible computational cost for any P(uj) as long as its value 
is not extremely large. Note that although these results 
are for perfectly conducting spheroids ctip = oo, calcula¬ 
tions for arbitrary finite Ctip can be done in the same way 
except one has to use Eq. (21) instead of Eq. (22). 

Let us now compare the obtained dependence of 
on Ztip with the limiting asymptotic behavior predicted 
in Sec. H. First, at Ztip <C a, the poles of the spheroid 
approach that of a sphere, as expected, see Fig. 4(a) and 
(c). The other limit is Ztip ^ L, where the point-dipole 
formulas (9a)-(9b) should apply. In Fig. 4 it is seen that 
the lowest eigenvalue of both shapes indeed have the cor¬ 
rect behavior. The intermediate regime a Ztip L 
is the most nontrivial one. We argue that in this regime 
function Pp{ztip) behaves as 

Pt) (■^tip) — oln(ztip/a) , Cl Ztip A , (38) 

with some coefficient c ~ 1 independent of L. To ar¬ 
rive at this formula we first find bounds on Pq using 
the following theorem. Consider two perfectly conduct¬ 
ing probes of different sizes. If the surface of one probe 
can be inscribed into the other, then the first probe must 
have a larger Pq. This statement is physically natural be¬ 
cause self-sustained oscillations around the smaller body 
require a larger compensation from the surface reflectiv¬ 
ity, cf. Sec. H. It can also be proven mathematically from 
the variational principle.To place bounds on Pq of 
the spheroid we can consider two other probes, a larger 
one and a smaller one. We get 

pcone,. < p. < pss,. ^ ( 39 ) 

where is the lowest pole of a cone with a vertex 

touching the sample and enveloping the spheroid; 
is the lowest pole of a spheroid of shorter length L = 
Ztip. It can be shown^® that /Jg™®’"'' ~ (I/tt) ln(ztip/a). 
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As for the smaller spheroid, the point-dipole formula 
should apply by order of magnitude, /Sq'*’"*' ^ 61n(ztip/a), 
cf. Eqs. (9a) and (35a). Since the functional form of 
these bounds coincides with Eq. (38) up to numerical co¬ 
efficients, we argue that /3f)‘(ztip) should obey the same 
equation as well. The graph shown in Fig. 4(a) is con¬ 
sistent with this prediction. However, due to numerical 
limitations L/ztip and Ztip/a could not be very large in 
our simulation and we could obtain only a crude estimate 
1 < c < 3 of the coefficient c. The poles /?jj of the in¬ 
plane polarizability, which are plotted in Fig. 4(c) as a 
function of a, also show crossovers among three regimes 
(short, long, and intermediate distances) and can be un¬ 
derstood in a similar way. 

The behavior of the residues R'^ is more difficult to 
analyze. At large distances Ztip ^ L, the residues of the 
spheroid approach those of the sphere [Eq. (10)]. At small 
distances, where the poles behave as ln/3^ ~ {2k + 3)a, 
the polarizability is determined by a large number ~ 1 /a 
of terms in the pole-residue series. The sum rule (4) im¬ 
plies that the sum of these dominant residues must be 
of the order of Xo for each shape. Indeed, the residues 
of the sphere, which have the form oc ka^a^ and 

oc (fc-|-l)(fc-|-2)a^a^ (Appendix B), obey this require¬ 
ment. The residues of the spheroid are always larger than 
those of the sphere, consistent with the higher Xo- The 
intermediate-distance behavior of Rk defies an obvious 
characterization. It is intriguing that at small distances 
only the residues are affected by the aspect ratio of the 
probe, while at large distances only the poles are altered. 

Information about the probe-sample coupling comple¬ 
mentary to the properties of the poles and residues can 
be obtained by examining the potential distribution of 
the polariton modes in real space. The examples for 
the T modes are depicted in Fig. 1(b). The potential 
is strongly peaked near the tip of the spheroid, demon¬ 
strating the localized nature of near-held coupling. Note 
that the number of times the potential changes sign along 
X is equal to k. 

IV. MOMENTUM-DEPENDENCE OF THE 
PROBE-SAMPLE COUPLING 

A simple physical picture of the s-SNOM that served 
as an important insight in the early days of the held 
and still remains popular today is the notion that the 
probe couples predominantly to momenta q 1/a. Ac¬ 
cordingly, the s-SNOM signal is collected from a very 
small region of size ~ a directly underneath the tip. 
Modern applications of s-SNOM to two-dimensional and 
layered systems require going beyond this oversimpli- 
hed picture because the g-dependence of the rehectiv- 
ity ?'p(<j-')<?) of such systems can be very sharp due to 
presence of dispersive collective modes (Sec. VII). Recall 
that for a momentum-independent rehectivity /3 [Eq. (2)], 
the poles and residues of the polariton eigenmodes are 
determined solely by the permittivity and geometry of 


the probe. Unfortunately, for a g-dependent rehectivity 
such a clean separation of the probe and sample prop¬ 
erties in the eigenproblem is not possible. While one 
can still dehne the eigenmodes by suitably modifying 
Eq. (3), the corresponding poles and residues will be, 
in general, complicated functionals of rp{ui, g). However, 
if the g-dependence of the rehectivity is weak, it can be 
treated as a perturbation, and the sample-independent 
resonant modes are retained. As we show in this Sec¬ 
tion, in this perturbative case one can precisely dehne 
the probe-sample coupling as a function of g and the 
‘dominant’ momentum as a function of Ztip- 

Consider a small g-dependent correction to the rehec¬ 
tivity: 

rp{q) =/3 + 6rp{q) . (40) 

The A:th pole /3/. of the probe-sample eigenmodes is a 
functional of rp. The key question is how this pole is 
affected by the nonlocal correction to rp. The answer can 
be written in terms of —G^(g, Ztip)) the hrst variational 
derivative of /3fe[rp(g)] with respect to rp: 

OO 

Sl^kiztip) = - j Gl{q, Ztip)5rp{q)dq . (41) 

0 

This is the desired relation to the leading order in Srp. 
A few general properties of function at g < 1/ztip can 
be established. First, this function decays exponentially 
at large g: 

G^(g,z«p)~e-2«^“p. (42) 

This is so because the probe-sample interaction is medi¬ 
ated by multiple rehections of evanescent waves (Sec. I) 
and the shortest distance such waves have to travel is 
2ztip. Next, it is easy to see that G/. is normalized: 

OO 

J G’/iq, ztip)dq = 1. (43) 

0 

Using a variation principle one can also show that for 
a perfectly conducting probe G^(g, Ztip) is nonnegative. 
Therefore, functions G^(g, Ztip) can be considered weight 
functions for the perturbation 5rp{q). To put it another 
way, this set of functions quantifies the momentum de¬ 
pendence of the probe-sample coupling. Below we show 
that the properties of these functions paint a much more 
nuanced physical picture than the naive idea that the 
coupling is maximized at a single momentum g ~ 1/a. 
However, if one insists on characterizing the entire dis¬ 
tribution of relevant momenta by a single number, the 
logical candidates are the average momenta 

OO 

€ = j Gl{q, ztip)qdq . (44) 

0 
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FIG. 5. (Color online) (a) The weight function Go"(g, atip) for the spheroid with L — 25a and etip = oo. The dashed line shows 
<7o“(-Ztip). (b) Go(q,Ztip) for several ztip, with circles indicating q(^. The solid lines are for the spheroid, the dashed lines are 
for the sphere. The spheroid is more sensitive to small q compared to the sphere, while both shapes are more sensitive to large 
q as Ztip decreases, (c) The first three G^{q) for Ztip = a, with solid circles indicating q^ . The number of nodes in G^{q) is 
equal to fc, while qj: increases with k. The logarithmic scaling of the horizontal axes is used to show the small-g structure more 
clearly. 


The idea is that unless G'^{q,Zt\p) has a complicated 
structure or a slow decay, q'^ should play the role of a 
characteristic momentum that determines fcth polariton 
pole Accordingly, we may expect that 1/q^ should 
give an improved estimate of the spatial resolution of the 
probe in the context of near-held imaging by s-SNOM. 
Interestingly, can be found by differentiating /3^(ztip): 

= (45) 


To obtain this formula consider hrst a sample with a q- 
independent rehectivity jS and let the probe-sample sep¬ 
aration be Ztip = z + dz. This system is equivalent to 
another one: the probe separated by ztip = ^ from a 
hctitious two-component medium composed of a vacuum 
layer of thickness dz plus the original sample. The sur¬ 
face reflectivity of such a two-component medium is q- 
dependent, rp(g) = so that it has the form (40) 

with Srp{q) = —2qj3dz. Evidently, such a 6rp{q) shifts 
the resonant pole from /3 = I3^{z) to /3 = l3^{z-\-5z), i.e., 
causes a differential change 5(3^ = {d(}'^/dz)dz. Substi¬ 
tuting these relations into Eq. 41, we get Eq. (45). Note 
that as rises more steeply with ztip for larger k, q^ 
increases with k. 

An equivalent description of the effect of a g-dependent 
perturbation is that it induces a correction to the surface 
reflectivity. The effective reflectivity is different for 
each k and v, 

OO 

= P-SP'j^ = j Gl{q, ztip)rp(g)dg . (46) 

0 

The corresponding polarizability x'^ is given by 


X 


= E 


Rl 


fc =0 


PI-PI'' 


eff 


(47) 


In the following we focus on function GQ{q, Ztip) because 
fc = 0 is the dominant resonance at all but very small 
Ztip. Actually, the large-distance limit of this function 
has the universal form 


Goiq, Ztip) zz 4zt%g^e , Ztip > L , (48) 

same for both ly. Equation (48) follows from Eqs. (11) 
and (41) and is consistent with the surmised large-g be¬ 
havior (42). As one can see, Eq. (48) gives GQ(g, ztip) 
that is normalized, nonnegative, and has a single max¬ 
imum at g = 1/ztip. The average momentum is gg ~ 
3/(2ztip). 

In the intermediate-distance regime functions 
G^(g, Ztip) are not expected to be universal. The 
specific example we treat in detail is again the con¬ 
ducting spheroidal probe. Combining Eq. (45) with the 
results of Secs. II and III, for the strongly elongated 
spheroid we obtain the following: 


1 

9o‘(2tip) 


(aztip)^/^ , 

2».,plog(iS) , 


3 


■^^tip ^ a, 

a < Ztip < L , 
Ztip > L. 


(49) 


Since the left-hand side has the physical meaning of the 
spatial resolution of the probe, we expect it to monoton- 
ically decrease as Ztip decreases. Therefore, the length 
scale L appearing on the second line of Eq. (49) should 
be of the order of L/31og(L/a). The presence of a 
large logarithmic factor log(ztip/a) in the intermediate- 
distance regime a ztip ^ L indicates that function 
Gj)‘(g, Ztip) has a considerable weight at g parametrically 
smaller than 1/ztip. In other words, a strongly elongated 
spheroidal probe senses electric fields beyond its imme¬ 
diate vicinity p < Ztip. (A similar point was made previ¬ 
ously in Ref. 28.) As L/a decreases, L comes close to a, 


























11 


and this intermediate regime disappears. For example, 
the sphere acts essentially as a local probe. 

The calculation of G'^{q,Ztip) for the spheroid can be 
done as follows. Applying the first-order perturbation 
theory to the linear system (30), one finds 


GkiQ, Ztip) 


ujH'ufc 

u|,Hufe 


(50) 


where H' is the matrix with elements 


For k > 0, G^{q, Ztip) has nodes as a function q at hxed 
Ztip. The number of nodes is equal to k, see Fig. 5(c). 
Apparently, at such q near-held coupling between oscil¬ 
latory charge distributions on the probe and the sample 
exactly vanishes. Therefore, small perturbations at such 
discrete q do not affect the kth resonant mode. Finally, 
although q^ increase with k for the reasons explained 
above, the maxima of G^ show the opposite trend, which 
is presently not understood. 


, Ott 

Hu, = -/,+ i(gA)4+i(gA)e-2«-- 
Q 


(51) 


V. FROM NEAR-FIELD POLARIZABILITIES 
TO FAR-FIELD OBSERVABLES 


Once the eigenvectors are found, e.g., as described in 
Sec. Ill, function G^(q,Ztip) can be readily computed. 

Our numerical investigation of G^(q, Ztip) was limited 
mainly to fc = 0 and i/ = J-. We observed that the eigen¬ 
vector components approximately followed the geomet¬ 
ric series (uo)j ~ tG The quotient t is somewhat larger 
than unity for small ztip. As Ztip increases, t becomes 
less than unity, so that the first component (uo)o dom¬ 
inates. Neglecting all other components and expressing 
the modified Bessel function Iz/ 2 {z) in terms of elemen¬ 
tary functions, we obtain the analytical approximation 
from Eqs. (50) and (51): 

G^((Z, Ztip) = ^ (gFcoshgF - sinhgF)' , (52) 

where cq is a normalization constant. At Ztip ^ L we can 
focus on the range of momenta less than 1/L because at 
larger q this function is already exponentially small. For 
such q the bracketed expression on the right-hand side 
can be replaced by {Fq)^/9 and Zp = ztip + L hy zup, 
which yields the asymptotic form (48). 

To examine small and intermediate distances we used 
the direct numerical evaluation of Uq and Go(g, Ztip)- As 
in Sec. Ill, we considered two aspect ratios: L/a = 25 
and Lja = 1. Only v = 1. part was studied. The re¬ 
sults for L/a = 25 are shown using the false color scale 
in Fig. 5(a). It can be seen that as Ztip decreases, both 
qi^(ztip) and the position of the maximum of Gj/(g, Ztip) 
as a function of q shift toward larger values. This implies 
that the probe becomes more sensitive to finer spatial 
features of the sample, as discussed above. The line plot 
of Gj/(( 3 ', Ztip) for several Ztip presented in Fig. 5(b) de¬ 
picts the same trend. The average momentum gj/ and 
the position of the Gp{q) maximum are of the same or¬ 
der of magnitude except at very short distances where 
q^ increases more rapidly as Ztip decreases. Note that 
Eq. (49) predicts that diverges at Ztip = 0. From 
Fig. 5(b) we also see that for the same Ztip the maximum 
of G[}-( 9 , Ztip) is found at q smaller by a factor of 3-10 
for the spheroid compared to the sphere. This confirms 
that the spheroid is much more sensitive to small in-plane 
momenta than the sphere, i.e., the response of a strongly 
elongated spheroid is affected by a relatively wide range 
of lengthscales. 


In order to apply our theory to simulation of s-SNOM 
experiments, we need to include a few more ingredients in 
our calculation. The hrst one is the so-called far-held fac¬ 
tor (FFF) F''(w). This factor accounts for the fact that 
the probe is illuminated not only by the incident wave 
but also by its rehection from the sample. In experiment 
P-polarized incident held is usually used, to take advan¬ 
tage of the high transverse polarizability of the probe. 
Assuming the sample surface is hat, uniform, and its lin¬ 
ear dimensions are much longer than the radian sphere 
diameter c/w, the rehection of the incident wave is de¬ 
scribed by the coefficient rp(gs,a;), where 

qs = — sin 9 (53) 

c 

is the in-plane photon momentum and 6 is the angle of 
incidence. Hence, the ratio of j^-component of the elec¬ 
tric held at the surface to that of the incident wave is 
I ± rp{qs,u!) for z/ = _L and ||, respectively. The FFF 
also takes into account that the held scattered by the 
probe reaches the detector in two waves: directly and 
after rehection from the sample surface. Usually, the 
backscattered held is measured. It has the in-plane mo¬ 
mentum —qs and therefore the same rehection coefficient 
rp{—qs,uj) = rp{qs,ui) as the incident wave. The total 
FFFs for this setup are given by 

F-'~{uj) = [I -b rp{qs, uj)Y sin^ 9 , (54a) 

Fll(tLi) = [I — rp((;s,a;)]^cos^0. (54b) 

The trigonometric factors on the right-hand side take 
care of conversion between the total electric held Lext of 
the waves and their _L, || components. Note that our as¬ 
sumption of the plane-wave illumination is not entirely 
realistic. In experiment, a focused Gaussian beam is typ¬ 
ically used, in which case the FFFs are effectively aver¬ 
aged out over a range of angles 9. Numerical apertures 
~ 0.4 are common. We must also stress that Eq. (54) 
should be modihed if the system studied by s-SNOM is 
nonuniform on scales shorter than c/lo. Typical examples 
include a small sample residing on some substrate^® or 
measurements done close to a boundary of two different 
materials. 
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Another point we have to discuss is signal demodu¬ 
lation. In the experiment the probe is made to oscil¬ 
late mechanically, which causes periodic variation of the 
probe-sample distance: 

= -Zo + Az (1 — cos (/?),(/? = nt. (55) 

The oscillation amplitude is typically Az = 20-90 nm, 
comparable to the radius of curvature a ^ 30 nm of the 
probe. The minimal approach distance zq > 0 can be 
equal to zero if the probe taps the sample. The tapping 
frequency is many orders of magnitude smaller than 
the laser frequency ui, and so the motion of the tip does 
not affect the electromagnetic response. Effectively, the 
experiment consists of measuring the scattered signal for 
many static conhgurations with different Ztip. The nth 
Fourier harmonic of the backscattered field is referred 
to as the demodulated signal s„. (Here we define as 
a complex number but in experimental literature it is 
common to discuss the amplitude and the phase of 
separately.) The primary purpose of demodulation is to 
suppress the far-held background signal created by rehec- 
tions from the body of the tip, the cantilever, etc. This 
background is large but depends on ztip very weakly (lin¬ 
early) and thus contributes predominantly to the n = 1 
harmonic. Unfortunately, demodulation strongly dimin¬ 
ishes the signal amplitude, making it more susceptible to 
experimental noise. In practice, n = 2 or 3 usually gives 
the best approximation of the true near-held signal. 

The demodulated signal is related to the polarizabili¬ 
ties x'^(a;,Ztip) we have been discussing in previous sec¬ 
tions by 


Sni^) = const X , (56) 

where x5((w) is the nth Fourier harmonic of 


TT 



0 


One more element of the experimental protocol is normal¬ 
ization. What is typically reported is s(((w) normalized 
against a certain reference material, e.g.. Si or Au: 

s-:((o;) = <(o;)/<'-f(..). (58) 

The normalization eliminates a number of physically un¬ 
interesting or poorly known factors, such as the constant 
in Eq. (56) that are related to the optical setup of the 
experiment. The FFFs may also be canceled if both the 
studied and the reference objects in the experiment are 
positioned nearby, so that the data for the two are taken 
at points no farther apart than the diameter c/w of the 
radian sphere. 

The last point we wish to draw attention to is that the 
absolute value of the minimum probe-sample distance zq 
[Eq. (55)] cannot be determined very accurately. There¬ 
fore, experimentalists have to measure the so-called ap¬ 
proach curve, which is the s-SNOM response as a function 


of Zq at a fixed frequency. They then identify the point 
Zq = 0 as a point where a qualitative change in behavior 
in S 2 or S 3 appears. The logic behind this procedure is 
that once the probe makes the mechanical contact with 
the sample, its oscillations become reduced in amplitude, 
marking an unambiguous change. A potential flaw of this 
argument is that sharp changes in s„’s may be generated 
by a rapid variation of electromagnetic coupling between 
the probe sample at short separation even before making 
mechanical contact. We will discuss this issue in more 
detail in Sec. VI. 


VI. CASE OF LOCAL REFLECTIVITY: 
ALUMINUM OXIDE 


In this and the following Sections we discuss the im¬ 
plications of our theory for near-held response of real 
materials. We choose bulk a-Al 203 , also known as sap¬ 
phire or corundum, as our hrst example of highly reso¬ 
nant material with a momentum-independent rehectiv- 
ity /3 [Fig. 2(a)]. Another material with these properties, 
silicon carbide, has been a subject of a recent s-SNOM 
study co-authored by two of the present authors.^® Mod¬ 
eling results based on the BEM showing good agreement 
with the data were also reported in that work. Realis¬ 
tic probe shapes and retardation effects have been taken 
into account in order to achieve that. The latter was nec¬ 
essary since the probe length 2L ^ 20 /rm in the experi¬ 
ments was in fact larger than the diameter of the radian 
sphere c/w ^ 11 /im. Here we do not aim for a perfect 
agreement with a particular experiment but instead wish 
to illustrate how the general theory of multiple eigen- 
modes formulated in the preceding Sections can generate 
novel features in near-held observables. We study mostly 
probes of an idealized spheroidal shape but examine some 
other shapes as well. We stay within the quasistatic ap¬ 
proximation but we will comment on retardation effects 
in Sec. VHI. 

We use the following momentum-independent model 
for the rehection coefficient of the uniaxial AI 2 O 3 crystal, 

/3(uj) = I , eeff(w) = y^, (59) 

Eeff + 1 


where Cp for p = o (ordinary) and e (extraordinary) axes 
is given by 





^jLO,p - 


*7jLO,pW 

*7iTO,pW 


(60) 


The optical constants of AI 2 O 3 reported in the litera- 
ture^*^di have slight variations, presumably because of 
different crystal purity and processing. In our calcula¬ 
tions we adopt the results of Ref. 41 at room temperature, 
reproduced in Table. H. (For simplicity, the weak oscilla¬ 
tor at wto,o = 634cm“^ is neglected.) Due to smallness 
of the optical phonon linewidths 7 p in this material, the 
near-held rehectivity of AI 2 O 3 can be as high as ~ 10. 
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FIG. 6. (Color online) Response of a perfectly conducting spheroidal probe with L = 25a and bulk AI 2 O 3 sample, (a) The 
false color plot of 2 tip)/a®. The bright curves correspond to the resonant modes, with fc = 0 mode having the lowest 

frequency, (b) The polarizability (absolute value and imaginary part) at Zt\p ~ 0.02a = 0.6 nm, the smallest distance in 
panel (a), (c) The absolute value of the demodulated polarizability lyal and scattering signal |s 3 | for the tapping amplitude 
Az = 50 nm and zo = 0.6 nm. The origin of the three peaks is discussed in the text, (d)-(f) The counterparts of panels 
(a)-(c) for the parallel component, = ||. The plots again reveal multiple resonances. However, the overall magnitude of 
the polarizability is greatly reduced, x'' ~ 10 “^x^i the resonances are more strongly bnnched near the surface phonon 
frequency ojsp = 818 cm“^. 


TABLE II. Parameters of the optical constant of a-ALOa 
used in our calculations, cf. Eq. (60). The frequency unit is 
1 cm“^. (From Ref. 41 for temperature T = 300 K.) 


p 

€00 

j 

UJjhO 

7 jLO 

UJjTO 

IjTO 

0 

3.05 

1 

908 

22.4 

569 

7.86 



2 

482 

2.96 

439 

3.23 



3 

387 

5.18 

384 

6.03 

e 

2.9 

1 

885 

21.6 

582 

4.17 



2 

481 

3.21 

482 

3.42 



3 

511 

1.42 

400 

4.68 


We start by studying the behavior of the probe polar¬ 
izabilities x'^ S'® ^ function of frequency w. In the mid- 
infrared range, the reflection coefficient /3 of AI 2 O 3 has 
a single peak centered at the surface-phonon frequency 
wsp = 818cm“^, depicted in Fig. 2(a). As uj approaches 
wsp from below, Re/3(a;) steeply rises. Equation (3) im¬ 
plies that whenever Re/3 is equal to a pole /3^, Imy*^ has 


a local maximum as long as the damping Im/3(a;) is not 
too large. The positions of three such underdamped res¬ 
onances are indicated schematically in Fig. 2(a). Thus, 
a single surface mode wsp of AI 2 O 3 may produce mul¬ 
tiple modes of the coupled probe-sample system. These 
localized eigenmodes (resonances) have been discussed at 
length in the preceding Sections. For example, they are 
depicted in Fig. 1(b) for the case of a spheroidal probe. 
Note that all the resonances are red-shifted from the fre¬ 
quency wsp. Since Ira. (3 increases as w approaches wsp, 
higher-order resonances are progressively more broad. 

The scenario above is described in terms of constant 
/3^. However, the poles are functions of Ztip, and so 
the frequency of each resonance shifts with ^tip- This 
is clearly seen in a false color plot of Imx-*-(a;, Ztip) 
[Fig. 6 (a)], where each mode creates a bright curve. All 
the curves are red-shifted from wsp but converge to it at 
large znp. The smallest ztip = 0.02a in Fig. 6 (a) is limited 
by the accuracy of our numerical calculation. Based on 
our analytical results we expect that at smaller ^tip the 
resonance curves are shaped as parabolas that approach 
WTO = 576cm“^ where Re/3 = 1, cf. Eqs. (5), (59), and 
(60). A horizontal line cut through Fig. 6 (a) taken at 
Ztip = 0.6 nm is plotted in Fig. 6 (b) along with the ab- 
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solute value of x'*"- The strongest peak in this plot cor¬ 
responds to the k = 0 mode. The multiple weaker peaks 
at higher frequencies are produced by fc > 0 modes. 

Next we consider the effects of demodulation on the 
s-SNOM signal, which can be understood as follows. As 
the probe oscillates, it spends most time at the mini¬ 
mum and maximum distances from the surface. One 
therefore expects peaks in Xn frequencies near those of 
X'^{zo,uj) and x "{zo + 2Az,uj). This gives two frequencies 
per each resonant mode. Actually, the number of observ¬ 
able peaks is smaller. Indeed, from Figs. 6 (a) and 6 (d) 
one can see that all the resonance curves modes should 
merge together at z = Zq + 2Az for typical Az ~ 50 nm. 
Hence, all the modes should produce a single common 
peak in the demodulated signal from such z. Further¬ 
more, while the peaks of x'^izo,^) are distinct, only a 
few strongest of them can survive the smearing effect 
of the demodulation. These expectations are supported 
by Fig. 6 (c), where we plot the normalized quantities 
X3(w,2o) = Xs/xT^ and S3 for ly = ±, assuming tapping 
amplitude Az = 50 nm, zq = 0.6 nm, and Si as the refer¬ 
ence material. In Fig. 6 (c) we see only three peaks. The 
peak at 650 cm“^ in jssj is produced by the dominant 
k = 0 mode. It has the same frequency as the k = 0 
peak in Fig. 6 (b). The second peak near 725cm“^ in 
I S 3 1 (which looks more like a shoulder in xa) is produced 
by the k = 1 mode at the Ztip = zq point. The re¬ 
maining third peak at 787 cm“^ is produced collectively 
by all the modes. A similar correspondence between the 
resonance curves of the polarizability function and the 
peaks in the demodulated signal is found in the = || 
component, cf. Figs. 6 (d)-(f). However, the lower k = 1 
peak is now very weak and is considerably blurred by 
the demodulation, Fig. 6 (f). Should we have considered 
a model with smaller dissipation, this and other high- 
order peaks would have been more clearly distinguishable 
in |s 3 |. Note that although the normalized and demodu¬ 
lated signal strength is comparable for the two v compo¬ 
nents, the polarizability for iz = || is orders of magnitude 
smaller so its contribution can be safely ignored. 

The discussion above pertain to horizontal cuts of 
X‘^(-Ztip, w). Taking a fixed-frequency (vertical) cut 
through Fig. 6 (a), and performing the demodulation for 
a range of minimum distances zg, one obtains the v = A 
approach curve for the scattering signal. An intriguing 
result of this analysis is the possibility of a nonmonotonic 
dependence of the approach curve on zq . The nonmono¬ 
tonicity is due to the crossing of the resonance curves of 
X‘*“ by the vertical line cut. Such crossings are found be¬ 
tween wto where Re/? = 1 and wsp where Re/? reaches 
its maximum. Near the low-frequency end of this in¬ 
terval, the k = 0 mode should be again dominant. It is 
expected to produce a peak in the approach curve, which 
would follow the same trajectory as the k = 0 curve in 
Fig. 6 (a), moving to larger Zq as w increases. Higher or¬ 
der modes should appear at frequencies closer to wsp and 
produce weaker peaks at smaller zq- The amalgamation 
of these peaks give rise to the nonmonotonicity of the 


approach curve. 

We show in Fig. 7(a) the S 3 approach curves for v = A 
for three frequencies. All the curves are normalized to 
their value at their left ends, zg = 0-6 nm. The approach 
curve for uj = 600 cm“^ decays monotonically with in¬ 
creasing zg because the cut at such ui does not cross any 
of the resonances. In the approach curve for 700 cm“^, 
a strong peak is seen at around 2 nm due to the cross¬ 
ing of the k = 0 resonance. The last approach curve, 
for 800 cm“^ contains a series of oscillations at small zg 
and a broad hump at large zg, due to the multiple reso¬ 
nance crossings. The approach curves for iz = |j plotted 
in Fig. 7(b) exhibit the same general trends as those for 
iy= A. 

The striking multi-peak spectra and anomalous non¬ 
monotonic approach curves we described above stem 
from the large rp of AI 2 O 3 and are not found in less 
resonant materials, see Sec. VHI and Ref. 19. This rich 
structure is also quite sensitive to the choice of Zg. If this 
parameter is too large, the peaks in the spectrum of the 
scattering signal merge together at w = wsp. If zg is too 
small, the resonance curves become very flat at w < wsp, 
so the corresponding peaks are smeared by demodulation 
and dwarfed by the wgp peak. Hence, there exists an op¬ 
timal value of Zg that allows one to resolve multiple peaks 
most clearly. For our AI 2 O 3 model this value is actually 
not too far from zg = 0.6 nm used in Fig. 6 . For example, 
the S 3 spectrum for a smaller zg = 0.06 nm is shown in 
Fig. 8 (dashed lines), where the k = 0 peak is much less 
pronounced while more higher order peaks become dis¬ 
tinguishable and form small steps. For even smaller zg 
the steps are further smoothed, eventually leaving only 
one peak near wsp- 

In addition to the value of zg, many other experimen¬ 
tal parameters and procedures can significantly alter the 
resultant spectrum. For instance, the experimental de¬ 
termination of Zg based solely on the s-SNOM approach 
curve can be inaccurate due to its possible nonmonotonic¬ 
ity, as discussed in the previous Section. It is generally 
incorrect to ascribe zg = 0 to the probe position at which 
the near-held signal has the highest amplitude. Such a 
protocol effectively yields a frequency-dependent zg. The 
difference from the spectra taken for a truly constant zg 
can be drastic, as illustrated in Fig. 8 . Conversely, the 
strong sensitivity of the near-held signal to the value of 
probe-sample distance may perhaps be used for a more 
accurate measurement of zg (although this may require 
knowing the curvature radius a and perhaps other details 
of the probe shape). 

The tapping amplitude Az is another parameter that 
affects the spectrum. When Az is small, the demodula¬ 
tion at nth order is roughly equivalent to taking the nth 
order derivative of x‘^(z:tip). Therefore, a material with a 
sharply varying approach curve yields a stronger demod¬ 
ulated signal than the material with a smoothly varying 
one. In our case the signal of AI 2 O 3 is normalized against 
Si, whose polarizability decays monotonically with Ztip 
[Fig. 9(b)]. As Az decreases, the polarizability of AI 2 O 3 
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FIG. 7. (Color online) Approach curves of IS 3 I for bulk AI 2 O 3 , normalized to the value at zg ~ 0.6 nm, for several characteristic 
frequencies. The lowest frequencies in both (a) and (b) are such that no resonance curves are crossed during the probe tapping 
motion. The approach curves are monotonic. For the middle pair of frequencies one crossing (of the fe = 0 resonance) does 
occur. At such crossing each approach curve has a peak. The last pair corresponds to the frequencies where IS 3 I is close to the 
maximum value in the spectral range studied. The approach curves have several peaks because of multiple resonance crossings. 



u; (cm ^) 


FIG. 8 . (Color online) Comparison of spectra of the quan¬ 
tity |s;^| using two different experimental protocols, at two 
different minimum approach distance zg. The value at each 
frequency is taken either from the maximum of the |s^| ap¬ 
proach curve (solid) or from a fixed zg (dashed). 


become increasingly oscillatory, while that of Si remains 
smooth. This results in the increased contrast of the de¬ 
modulated signal for the two materials for smaller Az 
[Fig. 9(a)]. 

Other than these controllable parameters, the scatter¬ 
ing signal is also dependent on the dielectric function of 
the probe itself. The calculation in the preceding discus¬ 
sion is done for a perfectly conducting probe, etip = 00 . 
In practice, near-held probes often have a Si core and a 
layer of metallic coating whose thickness ~ 20 nm can be 
smaller than the skin depth, i.e., the electric held pene¬ 
tration length of the metal. In this case, it may be more 
appropriate to set etip = esi ~ 11-7 in Eq. (21). Re¬ 
peating the calculations, we hnd that while qualitative 


features in the signal are retained, there are major quan¬ 
titative differences (Fig. 10). 

The discussion above shows that the rich structure of 
the s-SNOM signal found for the case of AI 2 O 3 sample 
is susceptible to many experimental parameters. (Retar¬ 
dation effects, discussed later in Sec. VIII, introduce fur¬ 
ther signihcant dependence on the probe geometry.) This 
presents a serious challenge to realistic modeling of s- 
SNOM experiments. On the other hand, these strong de¬ 
pendences arise only for highly crystalline material with 
low dissipation. For other, less resonant materials, the 
modeling can be quite robust, as discussed in Sec. VIII. 

VII. NONLOCAL REFLECTION FUNCTION 

The example material of the previous Section is a bulk 
crystal with a local (momentum independent) reflectiv¬ 
ity function. However, in many other systems studied 
through s-SNOM, including thin films, graphene, and 
multi-layered systems reflection is inherently nonlocal. 
Thus, it is imperative to study how the g-dependence of 
the reflectivity affects the probe-sample interaction. As 
mentioned in Sec. IV, a general description of such inter¬ 
action is challenging because the series representation of 
the polarizability 

(«) 

k 

has generalized eigenvalues Afc and residues TZk that are 
now complicated functionals of rp, cf. Eqs. (20) and 
(30). Still, we can attempt to analyze these expres¬ 
sions using the simple perturbation theory developed in 
Sec. IV, in which \k are computed from the poles of the 
g-independent theory, with corrections obtained by in¬ 
tegrating the weighting functions over the momentum. 
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FIG. 9. (Color online) (a) Spectra of |s 3 "| for zo = 0.6 nm and five different tapping amplitndes. The magnitnde of increases 
rapidly with decreasing Az. (b) The approach curves for AI 2 O 3 , Si 02 and Si, taken at frequencies corresponding to the 
largest peak in |s;^| (790 cm“^ for AI 2 O 3 , 1120 cni“^ for Si 02 , and an arbitrary uj for the frequency-independent case of Si). 
For AI 2 O 3 sample, |x^| shows multiple oscillations; for Si 02 sample, it has a single maximum at small Ztip', for Si, the approach 
curve decays monotonically with Ztip. As Az decreases, the approach curves become increasingly different. 




FIG. 10. (Color online) (a) x^ and (b) |x:^| and |s^| spectra of AI 2 O 3 for the case of a Si probe. All other parameters are the 
same as in Fig. 6 (b) and (c). The spectra retain the same structure as for a metallic probe (etip = 00 ). 


As shown below, this scheme produces qualitative agree¬ 
ment with the calculated s-SNOM response for graphene 
on bulk AI 2 O 3 . 


The AI 2 O 3 /graphene system has two collective modes 
(the upper and the lower one) that emerge from hy¬ 
bridization of the surface phonon of AI 2 O 3 , originally at 
CUSP « 750 cm“^ with the plasmon of graphene, uj{q) oc 
y/Jtvpq. (Coupling of substrate phonons to graphene 
plasmons has been probed by s-SNOM experiments with 
graphene/Si 02 systems. This and related work is re¬ 
viewed in Ref. 42.) The modes share the optical weight 
and exhibit a level-repulsion that causes both to be dis¬ 
persive. Both features depend on the chemical potential 
of graphene. Below we focus on the upper mixed mode 
and study its s-SNOM response for a range of fi, and com¬ 
pare the results with the perturbation theory method. To 
proceed, we need the formula for the reflectivity of the 
composite system. This formula is well-known (see, e.g.. 


Ref. 5) 


rp{q,U}) 


ey _ Co dTTo- 

kf fcg UJ 

ey eo dTrcr 

kf fcg UJ 


(62) 


Here ei = Ceff [Eq. (59)] is the permittivity of the lower 
half-space (AI 2 O 3 ), cq = 1 is that of the upper half-space 

(vacuum), — q^ is the z-component of the 

wave vector in medium j = 0 , 1 , and a = cf{q,uj+iT~^) is 
the conductivity of graphene, which we calculate within 

the random phase approximation^^d4 ^j^h a finite re¬ 
laxation time T~^ = 25cm“^. For q ^ uj/c, one finds 
~ iq and Eq. (62) reduces to 


Cl — 1 -I- 4jKq — 

rp[q,uj) = ^- - 

ei -I-1 -I- 47rg — 


(63) 
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FIG. 11. (Color online) (a) Collective mode dispersion of graphene/ AI2O3 system. The mode repulsion between the graphene 
plasmon and the AI2O3 surface phonons are evident. The false color stands for Re{[l + rp(g,a;)]g/feo}, which is a measure of 
power dissipation.^® This quantity is additionally raised to power 0.35 to reduce the contrast. The vertical dashed line marks 
q — 1/a. The faint curve just below uj = 500 cm“^ is a weak surface phonon'^'^ that we do not discuss. The chemical potential 
of graphene is /r = 1200 cm“^. (b) The solid curves are constant momentum q = 1/a line cuts through maps like (a) for several 
fi. The particular case of (a) is shown by the red curve (second solid curve from the right). The dashed curve is the same 
quantity computed for bulk AI2O3 without graphene, (c) and (d) 83 “(oj) computed using the g-dependent rp{q,u!) at 

Ztip = 0.02a and zo = 0.02a, respectively. Graphene chemical potentials fi for (b)-(d) are indicated in the legend of panel (c). 


which can be compared to Eq. (59). A convenient way 
to visualize the dispersion of the collective modes is to 
plot the imaginary part of rp(g, w), which represents the 
power dissipation in the system,^® as a false-color map. 
An example for fjL = 1200 cm“^ is shown in Fig. 11(a). In 
the low-g regime {Hvpq huj the lower bright 

curve is mainly the plasmon with dispersion lo oc y/fivpq, 
while the upper bright curve represents the dispersion¬ 
less AI 2 O 3 surface phonon. (The additional bright curve 
around w = 500 cm“^ is a weaker AI 2 O 3 surface phonon, 
which we do not discuss.) An increase in fj, leads to a 
steeper dispersion of the plasmon, which causes both 
hybrid modes to go up in frequency. Decreasing /j, 
has the opposite effect. Additionally, if drops below 
husp /‘2 ~ 380 cm“^, the upper mode falls into the in¬ 
terband transition region of graphene, which results in 
strong damping of the surface phonon. As we will see 
below, this causes the n = 300 cm“^ curve to look quali¬ 
tatively different from the rest in Fig. 11(b). Let us now 
discuss how the collective modes manifest themselves in 
the s-SNOM response. 

In the simplistic picture of the s-SNOM response, the 


probe-sample interaction is dominated by a single mo¬ 
mentum q = 1/a. If this assumption were accurate, we 
could set rp{q,u}) as /3{oj) and calculate the response us¬ 
ing the set of poles and residues established previously. 
We would then see peaks in the response generated by 
the upper hybrid mode. However, this crude approx¬ 
imation leads to higher peak frequencies than the cal¬ 
culation using the full rp(q,w), as seen in Figs. 11(b), 
11(c), and 11(d). Indeed, we have shown in Sec. IV 
that when the g-dependence in reflection is treated as 
a perturbation, each mode has its own range of sensi¬ 
tive momenta due to the inherent length scales in its 
potential distribution. The distributions change with an 
additional length scale — the tip-sample distance Ztip, 
so that the momentum weighting functions are depen¬ 
dent on Ztip as well, Gk = Gfc(g,Ztip). For each mode, 
these functions provide a means to average over momen¬ 
tum and find an effective g-independent sample reflec¬ 
tion /3|®(w), cf. Eq. (46), so that we can again apply the 
established pole-residue decomposition. Strictly speak¬ 
ing, the perturbative method cannot be applied here as 
the mixed mode may be strongly g-dependent. Even 
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so, we find a very reasonable agreement with the com¬ 
puted signal in the range of graphene chemical poten¬ 
tials /X = 600-1800 cm“^ that we study. We first consider 
peak frequencies in Im , which can be predicted by in¬ 
voking the resonance condition Re = (3k ■ For the 
lowest mode k = 0 and ztip = 0.02 a, there is a system¬ 
atic overestimate of the peak position by 20-30 cm“^ for 
/X = 600-1800 cm“^. The discrepancy is larger for higher 
/X at which the g-dependence of the upper hybrid mode 
is stronger. This discrepancy is due in part to the well- 
known general tendency of the first-order perturbation 
theories to overestimate the lowest eigenvalues. Next, for 
the k = 1 mode, the resonance condition is satisfied only 
for /X = 600cm“^ at w = 797 cm“^ and /x=1200 cm“^ at 
uj = 823 cm“^, which agree well with the smaller peaks 
in Im y-*-. At these frequencies Im /3f® are larger than the 
k = 0 case and the peaks have smaller magnitudes. For 
/X = 1800 cm“^, the resonance condition is not met and 
the very small peak at ox = 827 cm“^ in Im^-*- corre¬ 
sponds to where Re f3f^ is largest and thus closest to /3i. 
Finally, for fc > I, (3k is larger than Refor all frequen¬ 
cies and no peaks in Im x'^ are found. Seeing qualitative 
agreement in the polarizability, we proceed to analyzing 
the demodulated signal. 


As inferred in Sec. VI, the demodulated signal is 
strongest near the peaks in x-*-(zo,a;) and X'^i^o + 
2 A 2 ;,a;), where each peak is attributed to a resonant 
mode. For the dominant k = 0 mode, we find a set 
of corresponding peaks in ( 2 : 0 , 0 ;) at the same frequen¬ 
cies as those in w), as shown in Fig. 11(c) and 

Fig. 11(d). For the other set of peaks in the S 3 spectra we 
must consider how the situation is changed at zq 2Az. 
At such distances Ztip itself becomes the primary length 
scale and the sensitivity function Gk is shifted toward 
smaller momentum, where the upper mode has a flatter 
dispersion and its frequency is close to o;sp of the bulk 
AI 2 O 3 crystal. Therefore, this set of peaks should all ap¬ 
pear near o;sp, which is indeed the case. Repeating this 
procedure for the k = 1 mode, we find that the peaks 
it contributes are inseparable from the set of higher fre¬ 
quency peaks produced by the k = 0 mode as both have 
frequencies very close to wsp. Its contributions, however, 
alter the heights of these peaks. For instance, the k = 1 
peak is strongest in x-*-(zo,a;) for /x = 600cm“^ (among 
the four we used), so the high frequency peak in for 
this chemical potential has the largest relative magni¬ 
tude with respect to the low frequency peak. Thus we 
conclude the demodulated s-SNOM signal can be qual¬ 
itatively explained by the perturbative method, albeit 
with inaccuracy in the lower peak frequency. However, 
as we argued in Sec. VI, the lower frequency peak in the 
demodulated signal is mainly an artifact of the finite zq 
we are forced to use. If Zq were truly zero, only the peak 
near wsp would survive. 


VIII. MODEL-DEPENDENT EFFECTS 

The spheroid model differs from real s-SNOM probes 
in two important ways: i) the real probe resembles an 
inverted pyramid, ii) at infrared wavelengths, the length 
~ 10 /xm of the probe exceeds several times the diameter 
c/uj of the radian sphere. In previous literature it was as¬ 
sumed that these differences can all be neglected as the 
probe-sample interaction is focused around the apex of 
the probe [Fig. 1(b)], while contribution from the rest of 
the probe is canceled out during the process of demod¬ 
ulation and normalization. Hence, the exact shape of 
the probe is unimportant and the only relevant physical 
quantity is the apex radius of curvature a. Further, since 
the characteristic length scale a is well within the ra¬ 
dian sphere, a quasistatic description should suffice. This 
simplistic argument is backed by previous agreement be¬ 
tween the spheroid model and experiment. 28 

ever, we have shown that different probe shapes exhibit 
universal behavior only when Ztip/a is of the order of a 
few percent (cf. Fig. 4(a)). This range is much smaller 
than typical tapping amplitudes, so the majority of the 
s-SNOM response lies outside the universality regime and 
should indeed be probe shape dependent. Additionally, 
recent experiment and modeling have shown that a qua¬ 
sistatic formalism with ad hoc probe shapes is insufficient 
for highly resonant materials such as on silicon carbide.^® 

In this Section we re-examine these issues by exam¬ 
ining two materials, the highly resonant AI 2 O 3 and the 
dissipative Si 02 , and study the probe shape dependence 
of their response as well as electrodynamic corrections. 
We find that for dissipative materials shape dependence 
is weak and retardation effects are of less importance, 
so the spheroid model describes the s-SNOM experiment 
reasonably well. This explains the success of our model 
in reproducing the response of various materials in ex¬ 
periment. On the other hand, we find the response of 
resonant materials to be highly dependent on the probe 
shape and less well described within the quasistatic ap¬ 
proximation. For such materials a full electrodynamic 
treatment with the exact probe shape may be required. 
Common numerical methods suitable for electrodynamic 
treatment of light scattering by a spheroid near a surface 
include T-matrix method^®’^® and BEM.^® For the case 
of a sphere near a surface, the calculation of necessary 
matrix elements can be done efficiently using recursion 
technique similar to what we use here.^^ 

We consider the probe shape dependence and the retar¬ 
dation effects separately. To study the former, we simu¬ 
lated the s-SNOM signal of AI 2 O 3 samples obtained with 
spheroidal probes of different length. We also calculated 
(using BEM) the results for pear-shaped probes that 
may better mimic the inverted pyramids. As shown in 
Eig. 12(a), the signal for a pear-shaped probe is qualita¬ 
tively similar to that for the spheroid of the same length, 
but there are quantitative differences. For spheroids, we 
find that the signal strongly increases and the peak fre¬ 
quencies steadily decrease as the length of the probe in- 
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FIG. 12. (Color online) (a) The s-SNOM signal computed for AI 2 O 3 samples. The inset shows the probe shapes used 
(spheroidal, with L = 25a, and a pear-shaped). The two types of probes produce qualitatively similar but quantitatively 
different results, (b) Spheroids of longer length have drastically increased signal strength. The inset depicts the probe shape 
and the values of L/a used. Note that this quasistatic calculation neglects radiative damping and antenna resonances, see 
Sec. VIII. If included, such effects are expected to greatly reduce S;^. In all cases A 2 = 50 nm and a = 30 nm. The value of S 3 
is taken either from the maximum of the approach curves at each frequency (solid lines) or at the closest approach distance 
zq = 0 . 6 nm (dashed lines). 


creases at a fixed apex radius, as shown in Fig. 12(b). 
These features can be explained by the scale invariance 
of the problem. It implies that an increase in probe 
length is equivalent to a simultaneous decrease in tap¬ 
ping amplitude and the apex radius. The decrease in 
radius produces changes in both the poles and residues. 
The former explains the shift in peak frequencies. The 
latter is mostly canceled out by normalization. In turn, 
the decrease in tapping amplitude leads to a larger con¬ 
trast between the sample and the reference as discussed 
in Sec. VI [see Fig. 9(b)], so the signal strength is dra¬ 
matically increased. 

The strong probe-shape dependence found above seem 
to suggest that theoretical modeling of the s-SNOM ex¬ 
periments must always be done using the actual shape to 
be reliable. In fact, such a sensitivity to the probe shape 
pertains only to the highly-resonant, i.e., large /3 materi¬ 
als. In AI 2 O 3 this parameter reaches the maximum value 
of |,5| « 12, Fig. 2(a). For comparison, in Fig. 13(b) and 
Fig. 13(c), we show that the pear-shaped probe and the 
spheroid produced almost identical signals for amorphous 
Si 02 , a material with |/3| < 1.5. (For experimental stud¬ 
ies of this material see, e.g.. Refs. 6 and 28.) In this case, 
a factor of 16 increase in the probe length leads to only 
a doubled signal strength, compared to a nearly tenfold 
increase for AI 2 O 3 seen in Fig. 12(b). 

The results above are obtained within the quasistatic 
approximation. In reality, a probe half-length of 200o 
already exceeds the diameter c/lu of the radian sphere 
and one has to consider retardation effects. Naively, con¬ 
tributions from such effects should be eliminated by de¬ 
modulation, as they pertain to a length scale much larger 
than the tapping amplitude. However, we show that one 
contribution — the radiative damping — survives de¬ 


modulation. The radiative damping has an effect sim¬ 
ilar to a finite Im,d, i.e., the dissipation in the sample. 
Hence, for dissipative materials one can neglect radiative 
damping and still find reasonable agreement with exper¬ 
iment, while doing so for highly resonant materials may 
lead to qualitatively wrong results. Let us illustrate these 
statements using the simplest model for the probe — the 
point dipole. The electrodynamic interaction between 
the dipole and the sample with the dielectric constant e 
is given^° by a modified version of Eq. (12), 


g =c 


eA;g(g)-fcf(g) 

K(.g) efcg(g) + fcf(g) 


dq, 


(64) 


where the second fraction in the integrand is the full form 
of the reflectivity rp{q,uj). [It is obtained from Eq. (62) 
by setting a to zero.] 

Suppose u! and zup are fixed, then the above integral 
defines g'^ as a function of e, which is generally a complex 
number. Alternatively, g'^ is a function of /3 = (e—l)/(e-|- 
1). The integration domain Eq. (64) includes momenta 
q both inside and outside the light-cone. The radiative 
damping effect arises from the integration over former, 
i.e., the momenta q < k^. This part of the integral yields 
a negative imaginary contribution to g'^ ^ which shifts the 
pole of x'^ [Eq. ( 11 )] to the lower complex half-plane of 
13. The real parts of the poles also change but this is 
less conceptually important, see below. Consider now 
the remaining part of the integral, over momenta g > fcg. 
It is easy to see that if e = — g^/(g^ — w^/c^), then 


efco(g)-k fci(g) = 0, (65) 

so that there is a pole on the integration path. As a 
result, functions g'^ and x'^ have branch cuts at e G 
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FIG. 13. (Color online) (a) The reflection coefficient of Si 02 ‘*® has a larger imaginary part than AI 2 O 3 due to its inherent 
dissipation, leading to a weaker shape dependence in the s-SNOM signal, (b) The signal of the pear-shaped probe is very close 
to that produced by the spheroid, (c) Increasing the probe length leads to a much smaller increase in the signal strength. The 
overall shape of the spectrum is also preserved. All geometric parameters are the same as in Fig. 12. 


(- 00 ,- 1 ] in the complex e plane or equivalently at 
P S [1, 00 ) in the complex /? plane. These additional 
features are shown schematically in Fig. 2(c). The phys¬ 
ical origin of both the poles and the branch cut is quite 
clear. The discrete poles has been discussed at length 
in this article. They correspond to the polariton modes 
localized near the tip. Fig. 1(b). In turn, the branch cut 
corresponds to the continuum of delocalized surface po- 
laritons that exist without the probe. Indeed, Eq. (65), 
is the well-known equation for the spectrum of such ex¬ 
citations.^® 

Of the two features, the branch cut is not expected 
to affect the signal as the small-momentum contribution 
is greatly diminished by demodulation. Demodulation 
should also make less important the change in the real 
parts of the poles, because these real parts vary greatly 
with ztip on account of the tapping motion of the probe. 
However, the shift of the discrete poles away from the 
real axis is a qualitative change and its effects remain 
after demodulation. Our next objective is therefore to 
find this shift for the case of the spheroidal probe. 

A free standing spheroid has an effective polarizability 
given by 


AO.eff 


Xo 




( 66 ) 


to the lowest order in w/c when radiative correction is 
considered.Modifying A accordingly [cf. Eq. (33)], it 
is easily shown that this formula applies to our geometry 
as well. Namely, the s-SNOM polarizability corrected for 
the radiative damping is given by 


Xrad 




00 


E 



(67) 


Viewed in the complex P plane, this correction is equiva¬ 
lent to the shift of the poles Pk into the lower half-plane 


by —i(2/3)(a;/c)^i?fc (to the leading order in w/c). There¬ 
fore, both the radiative damping and the intrinsic dissi¬ 
pation in the sample play a similar role: they increase 
the distance from the poles to the curve traced by the 
surface reflectivity P of the sample as w varies [Eig. 2(c)]. 
Eor a dissipative material, the curve begins far from the 
poles, and so further increase in the distance produces 
little change. Conversely, for highly resonant materials 
the P{uj) curve passes close to the real axis, and so radia¬ 
tive damping may obscure or eliminate the fine features 
of the signals, such as multiple resonant peaks discussed 
in Sec. VI. It is worth noting however that while it may 
be important for s-SNOM in infrared or visible domains, 
the radiative damping should be rather weak in the (ex¬ 
perimentally more challenging) terahertz range, where 
typical s-SNOM probes would fit well inside the radian 
sphere. 

Finally, a class of retardation effects we have not ad¬ 
dressed here are antenna resonances arising when the 
length of the probe exceeds several times the diameter 
of the radian sphere. They give rise to additional peaks 
in the s-SNOM signal as a function of w. For most ma¬ 
terials such resonances are removed once the s-SNOM 
signal is normalized to a reference sample; however, for 
strongly resonant materials such as SiC and presumably 
also AI 2 O 3 we studied here, the cancellation is not com¬ 
plete.^® 


IX. DISCUSSION AND CONCLUSION 

Further progress in the s-SNOM and related areas 
of near-held microscopy requires a quantitatively reli¬ 
able procedure for determining the fundamental response 
function rp(( 3 ', ui) from the amplitude and phase of the s- 
SNOM scattering data, from which one can proceed to 
the next step of inferring the optical constants of the 
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studied sample. Typically, materials with a higher ab¬ 
solute value of rp{q,uj) produce a higher amplitude s- 
SNOM signal. However, the peaks in the s-SNOM signal 
are often red-shifted with respect to those in |rp(( 7 ,a;)| or 
Im rp{q,uj). 

Given additional information about the system, these 
inverse problems can be tackled by fitting the experi¬ 
mental data to the solution of the direct problem with 
a trial form of rp{q,uj) as the input.Unfortunately, 
the direct problem is also difficult to solve. The three- 
dimensional nature of this problem and the presence 
of widely different length scales make realistic simula- 
tions^^“®^ of s-SNOM experiments very computationally 
intensive. This led to popularity of simple ad hoc approx¬ 
imations known as the point-dipole^^“^® and the finite- 
dipole model,in which the actual charge distri¬ 
bution induced on the probe is approximated by a point¬ 
like image dipole or a combination thereof with additional 
point charges. 

The point-dipole model^ postulates that Eqs. (9a)-(12) 
that are rigorous in the asymptotic long-distance limit 
Ztip L remain qualitatively correct at much shorter 
Ztip if the input physical parameters are suitably renor¬ 
malized. Thus, the bare polarizabilities Xo become the 
adjustable parameters of the model. It is customary to 
assume that the in-plane polarizability yj] is negligible 
compared the out-of-plane one, which is taken to be 

Xo = ( 68 ) 

where a is of the order of the curvature radius of the 
tip. Another adjustable parameter®’®^ ^ ^ 1 specifies the 
position of the effective dipole inside the probe: 

Zp = ba + ztip . (69) 

Clearly, the point-dipole model accounts only for the 
sharp tip and ignores the body of the probe, as x" for the 
point-dipole in Eq. (11) is much smaller than Xo for a tip 
with L ^ a. If the point-dipole model were literally cor¬ 
rect, the radiating dipole of the probe in typical s-SNOM 
experiments would be so small that no measurable signal 
would be observed. 

The finite-dipole model improves upon the point-dipole 
one by including the missing antenna-like enhancement 
approximately. It assumes that the electric field of a 
spheroidal probe of length 2L is equivalent to that of 
several point charges of total zero charge that are po¬ 
sitioned inside the spheroid near both of its ends. Eor 
small ZtipiL, this model®^’^® yields the following func¬ 
tional form of the probe polarizability: 

pfdp / \ 

if j ■ 

where oc aL^. The finite-dipole model was shown 
to give a good qualitative agreement with s-SNOM data 
obtained for quartz, amorphous Si 02 , and SiC samples 
once parameters and are suitably adjusted.® 


Thus, the best fit to the data was achieved choosing the 
length 2L = 600 nm of the probe, which is about one 
third of the diameter cjuj ^ 1700 nm of the radian sphere. 
Interestingly, this is approximately the value of 2L in the 
quasistatic calculation for which one obtains, in the case 
of Si 02 sample, the same result for S 3 as one gets from the 
full electrodynamic calculation for a probe of a realistic 
(much longer) length.^® 

Agreement with the data notwithstanding, from the 
theory point of view Eq. (70) is unsatisfactory on at 
least three counts. First, does not follow the correct 
scaling L®/lnL as a function of L, thus underestimating 
the probe polarizability. Second, the constant term in 
Eq. (70) violates the general requirement that x —>■ 0 
as /3 — )■ 00, corresponding to the case when the applied 
field is screened completely by the induced charges in 
the sample. Third, goes to ~ 1.4 when z^p = 0. In¬ 
stead, all smooth probe shapes must behave as a sphere 
at Ztip ‘C a, and therefore yield /3o = 1 at Ztip = 0 . 
The fact that finite-dipole model violates these general 
requirements suggest its limited usability. Figure 14 is 
an illustration of how widely different the predictions of 
the four discussed s-SNOM models can be for the case 
of AI 2 O 3 . Additional examples of similarly large differ¬ 
ences for Si 02 and SiC samples can be found in previous 
works of the present authors and their collaborators.^®’^® 
All these examples compel us to conclude that the prior 
success of the point- and finite-dipole models in fitting ex¬ 
perimental data has to be due to insufficient range of the 
data, multitude of adjustable parameters, and also the 
demodulation and normalization procedures that mask 
the errors in both the functional form and the magni¬ 
tude of the calculated signal. 

Another way to explain the difference between the ear¬ 
lier ad hoc models and our GSM is as follows. For the 
case of a sample with a local reflectivity /3, the exact 
scattering problem of a dielectric probe near a surface 
reduces to a generalized eigenproblem,®^ that has an in¬ 
finite number of eigenmodes, as we discussed in Sec. I. In 
contrast, both the point- and the finite-dipole models at¬ 
tempt to approximate the infinite number of eigenmodes 
by a single one. 

Since the real-space potential distribution of the eigen¬ 
modes [Fig. 1(b)] depends on the shape and size of the 
probe and probe-sample distance but not on fi, we can 
describe interaction of the probe with an arbitrary sam¬ 
ple efficiently using the precalculated basis of such eigen¬ 
modes. This allows one to use our GSM approach to 
model s-SNOM response for a wide range of materials. 
However, calculations for realistic probe shapes are not 
always practical. In search of a broadly applicable yet 
simple model, we have chosen the prolate spheroid to be 
our probe shape, as it captures the essential features of 
the actual probes — a sharp apex and a strongly elon¬ 
gated shaft. We quantified the eigenmodes of the probe- 
sample system in the form of poles and residues of the 
polarizability functions x^ (Table I), allowing an expe¬ 
dient, in fact, instantaneous calculation of the s-SNOM 
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FIG. 14. (Color online) The spectrum of the probe polariz¬ 
ability Ix^l for AI2O3 sample according to four different mod¬ 
els. The point- and the finite-dipole models each predict a sin¬ 
gle peak in |x^|- The calculations for spherical and spheroidal 
probes reveal multiple peaks. The sphere and the point-dipole 
models produce ~ 1- The L = 25a spheroidal probe 

yields x^ ~ 10^-10"^; the finite-dipole of the same L gives 
about an order of magnitude lower x^- These dramatic dif¬ 
ferences in both the form and the absolute magnitude of x^ 
can however be significantly reduced in the usually reported 
S 3 ", the normalized demodulated signal. 


reponse. The point-dipole, finite-dipole and other ad hoc 
models no longer have the advantage of computational 
speed and should now be considered obsolete. 

Recent work^® has shown that in the strong-coupling 
regime of the probe-sample interaction a fully electro¬ 
dynamic treatment using the BEM and realistic probe 
shape is necessary in order to reproduce the measure¬ 
ments. This regime is realized experimentally®’^® when 
using samples of SiC, a material for which |/3| can be as 
high as 15. The same considerations apply for AI 2 O 3 
for which |/3| can reach 12, see Fig. 2(a). Our GSM 
theory gives analytical insight into near-held response of 
such materials. We have shown that due to simultaneous 
excitation of multiple eigenmodes, novel features of the 
s-SNOM signal such as multi-peaked spectra and non¬ 
monotonic approach curves can appear. These features 
are however very sensitive to experimental parameters 
such as tapping amplitude, minimum approach distance, 
and even the data collection protocol. Retardation ef¬ 
fects, especially radiative damping can also qualitatively 
alter the signal and must be considered. In order to 
observe the predicted anomalous approach curves and 
multi-peak spectra, it may be necessary to make efforts to 
minimize the radiative damping, which requires working 
with shorter probes or at lower frequencies. In contrast, 
in the weak- and moderate-coupling regimes, which are 
relevant for the vast majority of samples, the lowest-order 
eigenmode is dominant. Hence, the approach curves 
should be monotonic in ztip, while the spectra should 
be mostly insensitive to experimental details and retar¬ 
dation effects. This is the regime where our spheroidal 


probe model can be used with the greatest confidence. 

Our GSM theory also applies to a more complicated 
problem where the sample reflectivity is nonlocal, i.e., 
momentum-dependent. Here the salient advantages of 
our method are two-fold. First, in the case of a weak 
nonlocality, our GSM provides a mapping of the nonlocal 
problem to a local one. Thereby the sample-independent 
eigenmode decomposition is retained, providing an in¬ 
tuitive interpretation of the scattering signal. Second, 
our numerical algorithm (see Supplementary online ma¬ 
terials) is much more efficient than the standard BEM 
because the number of necessary matrix element calcu¬ 
lations scales linearly instead of quadratically with the 
matrix size. It will be worthwhile to compare the ac¬ 
tual computational speed of our algorithm with that of a 
recently developed and significantly more efficient BEM 
that utilizes pre-calculated matrix elements.^® 

We hope that the improved physical understanding of 
near-held probe-sample coupling enabled by the gener¬ 
alized spectral method advanced in this work as well as 
the numerical procedures we developed for its implemen¬ 
tation can be of use for modeling and analysis of future 
s-SNOM and other near-held experiments. 
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Appendix A: The electrostatic problem of a 
spheroidal probe 

The electric held created by a spheroidal object is most 
conveniently described in the prolate spheroidal coordi¬ 
nates (^, rj, (j)) where the origin of the coordinate system is 
located at the center of the probe, as shown in Fig. 3. The 
relationships to the cylindrical polar coordinates (p, (j), z) 
are 

Z = F ^Tj , p = _ 7 y 2 _ 

In the spheroidal coordinate each spatial position is spec- 
ihed by ^ G [l,c»), rj G [— 1 , 1 ], and (j) is the usual 
azimuthal angle. Contours of constant ^ are a series 
of concentric spheroids centered at the origin, with the 
major axis along the z direction and common foci at 
z = ±F. For each such spheroid, ^ is equal to the ra¬ 
tio of its major semi-axis and focal length. We consider 
the case when the surface of the probe coincides with 
one of the spheroidal surfaces ^ = ^0 = L/F, where L 
is the half-length or major semi-axis of the probe. Re¬ 
lated quantities such as the minor semi-radius W of the 
probe or the radius of curvature a at the apex are given 
hy W = and a = (L. 
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It is well known that Laplace’s equation = 0 has 
separable solutions in the prolate spheroidal coordinates. 
In particular, we are interested in solutions outside a 
spheroidal probe that decay at large Their most gen¬ 
eral form is written in terms of a linear combination of 
spheroidal harmonics as follows: 

(A2) 

where are coefficients to be determined from bound¬ 
ary condition, P["(C) is the associated Legendre polyno¬ 
mial defined on the interval [—1,1] and Pp^iC) and Q™(C) 
are the associated Legendre function of the first kind and 
second kind (See, e.g. Ref. 59) with 

= min(^, ^o), ?> = max(^, ^q) • (A3) 

With the above definition of Eq. (A2) covers both 
inside and outside the surface of the probe at 

For the geometry considered in this paper [Fig. 1(a)], 


the total potential $ can be written as 

^ — ^0 ffi ^plane ^sphd 5 (A4) 

where 

$o(r) = -Eor (A5) 

is the potential of the external uniform field, and $piane 
is the potential due to charges in the sample, which can 
be decomposed into evanescent plane wave as 

<i>piane(r) = J B {q)e-'^^ ^ , (A6) 

where the position vector r = p + zz is broken up into its 
cylindrical polar coordinate components. We determine 
B{q), as well as R™ from boundary conditions. 

To do that we quote two well-known mathematical re¬ 
sults: the decompositions of evanescent plane waves in 
terms of spheroidal harmonics and vice versa. The first 
reads 




l—O m——l 


with (j) and (pq being the azimuthal angles of p and q, respectively. The reverse is: 


gr(c)pr(^) = ^ 


{I — my. J q \ qF 




y2 

iqp+qz^-im{rj>-rj>g) ^ 9 

47r^ ’ 


(A7) 


(AS) 


where Ii^{z) is the modified Bessel function of the first 
kind.^^ These two relations follow easily from addition 
theorems of general Legendre functions such as those in 
Ref. 59 and 60. 

Near the sample surface z = —Zp, the boundary con¬ 
dition is 

$ = R(q)e"«^p-f4phd(q,-2:p) oc -rp(g)e"'?^p , (A9) 
where we use the notation 


summarized by first dehning an infinite matrix: 

OO 

Hu'=2n J rp(g)/p+i(<7F)/,+ i(qF)e-2«--^, (A13) 

0 

whose elements are integrals of rp{q). Then for each 
integer m a quantity related to R™, a column vector 
defined by: 


f{Q,z)= J /(^)e (AlO) 

for a partial Fourier transformation. Eq. (A9) implies: 

B{q) = -rpe"'J"’p$sphd(q, ^p) ■ (All) 

The other boundary condition for a uniform spheroidal 
probe with dielectric constant etip is: 




{I — to)! F 


(A14) 


is the solution to the linear system of equations 


^ (A"* - H)^^, A^i, = b^i. (A15) 




9$ 




(A12) 


Boundary conditions in Eqs. (All) and (A12) with the 
decompositions in Eqs. (A7) and (A8) allow one to com¬ 
pute the unknown coefficients RJ". The result can be 


The diagonal matrix elements A™;// = are defined 

by Eq. (21), and the numbers on the right-hand side of 
the equations are given by 


4(1 +to)! 

° I = F71- 

3 (1 — m)\ 


(A16) 
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where = —E^ and = {E^ T iEy)/2. Thus, the 
form of matrix A is determined completely by the ge¬ 
ometry of the probe (in terms of ^o) and its dielectric 
constant Ctip, while H describes the interaction between 
the sample reflection function rp and the momentum se¬ 
lectivity of the modes. The column vector b describes 
the uniform external field. 

The coefficients can be obtained directly from 
Eq. (A14) after one solves for A^i from Eq. (A15). But 
for the purpose of determining the induced probe dipole 
moment, only |to| < 1 cases are important. By examin¬ 
ing the asymptotic behavior of the electrostatic potential 
$(r), one obtains the total dipole moment of the spheroid 
probe. Its Cartesian components are related to the com¬ 
ponents of (A™); by 


Psp,0 — Psp, 2 : 


p3 

^ ’ ^sp,l 


Psp,x Wsp,y _ 

2 ^ “ T ^ 

(A17) 


Appendix B: The spherical probe limit 


The spheroidal probe model presented in Appendix A 
is quite general and can be a good model for tips of any 
aspect ratio L/a. Here we explore a particular limit of 
F —0 and ^ oo while keeping the product —>• a 

constant. This corresponds to the problem of a spherical 
probe of radius a. The derivation in Appendix A sim¬ 
plifies to that in Sec. 4.1 of Ref. 20, and by using the 
following asymptotic forms of various special functions: 



QTiO 

pno 


22 '+U! [qEY 

{21 + 1)1 V 2 y ’ 

(Bla) 

2Hl{l + m)l 
^ ’ (21 + 1)1 ^ 

(Bib) 

( 20 ! w 

2H\{l-m)r ’ 

(Blc) 


one can show that the decompositions Eqs. (A7) and (A 8 ) 
reduce to Eq. (4.9) and (4.10) of Ref. 20. The character¬ 
istic equation Eq. (A15) for B”^i reduces to: 


yf Sg, {l + iy. ^ 

{I + - m)\ 


B'^ 


(B2) 


where 

= (_)'+-(!+ m)!^^^^F'+iF™z (B3) 

is similarly related to the induced charge distribution of 
the probe. 


ai 


;(etip - 1) 

^Ctip + 1) -f 1 


(B4) 


is the multipole polarizability of the probe, and 

OO 

Pl = lJ rp{q)q^e-^'^‘^»dq (B5) 

0 

is the integral that characterizes the interaction between 
the spherical probe and the sample with do = a + zup. 
Eq. (B2), the characteristic equation for a spherical 
probe, is derived in Ref. 20 as Eq. (4.20). The solution to 
Eq. (B2) has some of the same properties as the spheroid 
case: B'^i = 0 for all I and |m| > 1 ; B^ is related to the 
charge distribution due to the z component of the elec¬ 
tric field and B^^ are related to the charge distribution 
due to the x-y component of the electric field. 

For the case of q independent rp{q, ui) = /3(a;), in which 
the integrals F; reduces to: 


(2do)'+i 


(B 6 ) 


there is an exact solution to the spherical characteristic 
equation. Let 


= arccosh — = arccosh fl + ^ (gj) 

a \ a / 


be a dimensionless parameter that characterizes the 
sphere-to-sample distance relative to its size, and let 


O’fc 


(/3;a) = ^ 


m—0 

Using the following quantities: 
Po = XoEz , Xo = a^ , 

90 

a \ 

/ sinh a 


( 2 m-fl)^ 

g(2m+l)a _ a 


= — ( cosh a — sinh a — ) , 
a \ Op) 


Pn = Po/3" 


qn = 


\sinh(n -I- l)a 
/3" sinh a 


sinh(n + 1)q: 


9o - 


Po sinh na 
a sinh(n + l)a 


it can be shown that for a metallic sphere 

OO 

Bi = {-Y X! ^"(^0 “ “ Pnl{do - dn) 

where 


^-1 


n—O 


(B 8 ) 

(B9a) 

(B9b) 

(B9c) 

(B9d) 

(BIO) 


, , sinh na 

do dji = u . , . \ 

sinh(n + 1 )q: 


(Bll) 


The physical meaning of Eq. (BIO) becomes clear when 
one treats the problem with method of images (Fig. 15). 
Suppose that an external electric field F^ would have 
induced a bare dipole moment pg in the sphere. This 
would induce an image dipole in the dielectric half-space. 
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FIG. 15. The method-of-images solution of the problem of a 
metallic sphere above a dielectric half-space with the external 
field normal to the interface. The method involves an infinite 
series of dipoles and point charges located inside the 
sphere at distances above the interface. The total charge 
inside the sphere is zero. 


with (Tfe given by Eq. (B8). 

The above analysis resulting in Eq. (B12) is for the case 
where the electric field is perpendicular to the sample. 
Eor the case where the electric field is parallel to the 
sample, the analysis is simpler in that the positions of 
the image dipoles and their strength are the same, but 
no image point charges are present. Therefore, in this 
polarization: 


^sph,|| 

Xo 


total 

^xy . 1 3 

—-— = Sinll ( 
Po 


: (0-2 - (To) ■ 


(B13) 


Both Eq. (B12) and (B13) conforms to our earlier as¬ 
sertion that X has the form of Eq. (3): 


X' 


sph.ll _ 


osph,|| 

k=0 Pk 


,sph,|| 


-/3 


X' 


oo psph,_L 

sph,_L _ _ 

~ 2-^ osph,_L _ o ■ 

fc=0 Pk P 

(B14) 


Eor horizontal electric fields, is singular whenever 

(72 or (To is, so that 


osph.ll _ J2k+3)a 

Pk — ^ 


(B15) 


which would in turn induce an image dipole and an im¬ 
age charge in the sphere. The position and strength of 
each successive image dipole and charge can be solved 
by recursion. Setting the sample surface to z = 0, the 
position of the center of the sphere is at 2 ; = (io- At each 
position z = dn given by Eq. (Bll) there is a point dipole 
Pn and a point charge . Charge qo is determined by the 
neutrality condition ^^qn = Oj which yields Eq. (B9b). 
Summing up all the contributions to the total dipole mo¬ 
ment from both the dipoles and the point charges inside 
the sphere, we get: 


The corresponding residues are: 

^^ph.ll =4(yt + i)(/s + 2)xosinh3a. (B16) 

Eor electric fields perpendicular to the sample, the paren¬ 
thesis in Eq. (B12) vanishes at each so they are 

not poles of x®p''A. Instead, occur at the zeros of 

(To which has no simple analytic form. The poles of the 
two polarizations, however, interleave: 


^sph.J_ 

Xo 


P 


total 

z 


Po 


1 

Po 


00 

[pn + qnidn 

n—0 


= 2 sinh^ a 




(B12) 


osph.ll ^ osph,_L ^ osph.ll 

Pk ^ Pk ^ Pk+1 ■ 


(B17) 


1 H. Wheeler, Proc. IRE 47, 1325 (1959). 

^ F. Keilmann and R. Hillenbrand, Phil. Trans. Roy. Soc. 

London, Ser. A 362 , 787 (2004). 

® J. M. Atkin, S. Berweger, A. C. Jones, and M. B. Raschke, 
Adv. Phys. 61 , 745 (2012). 

^ S. Amarie, T. Ganz, and F. Keilmann, Opt. Exp. 17 , 
21794 (2009). 

® Z. Fei, G. O. Andreev, W. Bao, L. M. Zhang, A. S. 
McLeod, G. Wang, M. K. Stewart, Z. Zhao, G. Dominguez, 
M. Thiemens, M. M. Fogler, M. J. Tauber, A. H. Cas¬ 
tro Neto, C. N. Lau, F. Keilmann, and D. N. Basov, Nano 
Lett. 11 , 4701 (2011). 

® S. Amarie and F. Keilmann, Phys. Rev. B 83, 045404 


( 2011 ). 

^ F. Huth, M. Schnell, J. Wittborn, N. Ocelic, and R. Hil¬ 
lenbrand, Nat. Mat. 10, 352 (2011). 

® H. Ma and J. Levy, Nano Lett. 6, 341 (2006). 

® M. M. Qazilbash, M. Brehm, B.-G. Chae, P.-C. Ho, G. O. 
Andreev, B.-J. Kim, S. J. Yun, A. V. Balatsky, M. B. 
Maple, F. Keilmann, H.-T. Kim, and D. N. Basov, Science 
318 , 1750 (2007). 

H. Zhan, V. Astley, M. Hvasta, J. A. Deibel, D. M. Mittle- 
man, and Y.-S. Lim, App. Phys. Lett. 91 , 162110 (2007). 
M. M. Qazilbash, M. Brehm, G. O. Andreev, A. Frenzel, 
P.-C. Ho, B.-G. Chae, B.-J. Kim, S. J. Yun, H.-T. Kim, 
A. V. Balatsky, O. G. Shpyrko, M. B. Maple, F. Keilmann, 




















26 


and D. N. Basov, Phys. Rev. B 79, 075107 (2009). 

A. Frenzel, M. M. Qazilbash, M. Brehm, B.-G. Chae, B.-J. 
Kim, H.-T. Kim, A. V. Balatsky, F. Keilmann, and D. N. 
Basov, Phys. Rev. B 80, 115115 (2009). 

A. C. Jones, S. Berweger, J. Wei, D. Cobden, and M. B. 
Raschke, Nano Lett. 10, 1574 (2010). 

K. Lai, M. Nakamura, W. Kundhikanjana, M. Kawasaki, 

Y. Tokura, M. A. Kelly, and Z.-X. Shen, Science 329, 190 

( 2010 ). 

M. P. Nikiforov, S. C. Kehr, T.-H. Park, P. Milde, U. Zer- 
weck, C. Loppacher, L. M. Eng, M. J. Therien, N. Engheta, 
and D. Bonnell, J. Appl. Phys. 106, 114307 (2009). 

Z. Fei, A. S. Rodin, G. O. Andreev, W. Bao, A. S. 
McLeod, M. Wagner, L. M. Zhang, Z. Zhao, M. Thiemens, 
G. Dominguez, and et ah. Nature 487, 82 (2012). 

J. Chen, M. Badioli, P. Alonso-Gonzalez, S. Thongrat- 
tanasiri, F. Huth, J. Osmond, M. Spasenovic, A. Centeno, 
A. Pesquera, P. Godignon, A. Zurutuza Elorza, N. Ca¬ 
mara, F. J. G. de Abajo, R. Hillenbrand, and F. H. L. 
Koppens, Nature 487, 77 (2012). 

S. Dai, Z. Fei, Q. Ma, A. S. Rodin, M. Wagner, A. S. 
McLeod, M. K. Liu, W. Gannett, W. Regan, K. Watanabe, 

T. Taniguchi, M. Thiemens, G. Dominguez, A. H. C. Neto, 
A. Zettl, F. Keilmann, P. Jarillo-Herrero, M. M. Fogler, 
and D. N. Basov, Science 343, 1125 (2014). 

A. S. McLeod, P. Kelly, M. D. Goldflam, Z. Gainsforth, 
A. J. Westphal, G. Dominguez, M. H. Thiemens, M. M. 
Fogler, and D. N. Basov, Phys. Rev. B 90, 085136 (2014). 
G. Ford and W. Weber, Physics Reports 113, 195 (1984). 

M. Agranovich, B. Katsenelenbaum, A. Sivov, and 

N. Voitovich, Generalized method of eigenoseillations in 
diffraction theory (Wiley-VCH, Berlin, 1999). 

D. J. Bergman, Phys. Rep. 43, 377 (1978). 

D. J. Bergman, Phys. Rev. B 19, 2359 (1979). 

R. Hillenbrand and F. Keilmann, Phys. Rev. Lett. 85, 3029 

( 2000 ). 

T. Taubner, F. Keilmann, and R. Hillenbrand, Nano Lett. 
4, 1669 (2004). 

J. Aizpurua, T. Taubner, F. J. G. de Abajo, M. Brehm, 
and R. Hillenbrand, Opt. Express 16, 1529 (2008). 

N. Ocelic, Quantitative Near-Field Phonon-Polariton Spec¬ 
troscopy, edited by G. Abstreiter, M.-C. Amann, M. Stutz- 
mann, and P. Vogl, Selected Topics of Semiconductor 
Physics and Technology, Vol. 89 (Walter Schottky Insti- 
tuts der Technischen Universitat Miinchen, 2007). 

L. M. Zhang, G. O. Andreev, Z. Fei, A. S. McLeod, 
G. Dominguez, M. Thiemens, A. H. Castro-Neto, D. N. 
Basov, and M. M. Fogler, Phys. Rev. B 85, 075419 (2012). 
R. W. Rendell and D. J. Scalapino, Phys. Rev. B 24, 3276 
(1981). 

P. K. Aravind and H. Metiu, J. Phys. Chem. 86, 5076 
(1982). 

P. K. Aravind and H. Metiu, Surf. Sci. 124, 506 (1983). 
P. Bobbert and J. Vlieger, Physica A 147, 115 (1987). 

J. D. Caldwell, A. V. Kretinin, Y. Chen, V. Giannini, 

M. M. Fogler, Y. Francescato, C. T. Ellis, J. G. Tis- 
chler, C. R. Woods, A. J. Giles, M. Hong, K. Watanabe, 


T. Taniguchi, S. A. Maier, and K. S. Novoselov, Nat. 
Comm. 5, 5221 (2014). 

Z. Sun, A. Gutierrez-Rubio, D. N. Basov, and M. M. 
Fogler, arXiv: 1502.06524 (2015). 

C. F. Bohren and D. R. Huffman, Absorption and Scatter¬ 
ing of Light by Small Particles (Wiley, New York, 2004). 

®® MATLAB Release 2012b, The MathWorks, Inc., Natick, 
Massachusetts. 

P. M. Morse and H. Feshbach, Methods Of Theoretical 
Physics (Feshbach Publishing, Minneapolis, 1981). 

M. M. Fogler, “Electrostatic polarizability of a conical 
metallic particle in a contact with a dielectric half-space,” 
(2013), unpublished. 

®® F. Keilmann and R. Hillenbrand, “Nano-optics and Near- 
held Optical Microscopy,” (Artech House, Norwood, 2009) 
Chap. 11: Near-Field Nanoscopy by Elastic Light Scat¬ 
tering from a Tip, pp. 235-266, edited by A. Zayats and 

D. Richards. 

A. K. Harman, S. Ninomiya, and S. Adachi, J. Appl. Phys. 
76, 8032 (1994). 

Zeidler, S., Posch, Th., and Mutschke, H., Astron. Astro- 
phys. 553, A81 (2013). 

D. N. Basov, M. M. Fogler, A. Lanzara, F. Wang, and 
Y. Zhang, Rev. Mod. Phys. 86, 959 (2014). 

B. Wunsch, T. Stauber, F. Sols, and F. Guinea, New J. 
Phys. 8, 318 (2006). 

E. H. Hwang and S. Das Sarma, Phys. Rev. B 75, 205418 
(2007). 

T. Wriedt and A. Doicu, Optics Gommun. 152, 376 
(1998). 

'^® A. Doicu, Y. Eremin, and T. Wriedt, Optics Gommun. 
159, 266 (1999). 

P. Bobbert and J. Vlieger, Physica A 137, 209 (1986). 

A. Kuchkova and K. Navratil, Appl. Spectrosc. 48, 113 
(1994). 

'^® G. Ford and W. Weber, Surface Science 109, 451 (1981). 
®® A. Wokaun, J. P. Gordon, and P. F. Liao, Phys. Rev. Lett. 
48, 957 (1982). 

A. Moroz, J. Opt. Soc. Am. B 26, 517 (2009). 

J. A. Porto, P. Johansson, S. P. Apell, and T. Lopez-Rfos, 
Phys. Rev. B 67, 085409 (2003). 

R. Esteban, R. Vogelgesang, and K. Kern, Nanotechnology 
17, 475 (2006). 

R. Esteban, R. Vogelgesang, and K. Kern, Opt. Express 
17, 2518 (2009). 

S. Amarie and F. Keilmann, Phys. Rev. B 84, 199904(E) 

( 2011 ). 

®® B. Hauer, A. P. Engelhard!, and T. Taubner, Opt. Express 
20, 13173 (2012). 

J. Renger, S. Grafstrom, L. M. Eng, and R. Hillenbrand, 
Phys. Rev. B 71, 075410 (2005). 

A. Cvitkovic, N. Ocelic, and R. Hillenbrand, Opt. Express 
15, 8550 (2007). 

®® 1. S. Gradshteyn and 1. M. Ryzhik, Table of Integrals, 
Series, and Products, 6th ed., edited by A. Jeffrey and 
D. Zwillinger (Academic, San Diego, 2000). 

®® E. W. Hobson, Proc. London Math. Soc. s2-29, 355 (1929). 


